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Abstract 

Methods for fusing two computer vision methods are discussed and several example 
algorithms are presented to illustrate the variational method of fusing algorithms. The 
example algorithms solve the photo-topography problem; that is, the algorithms seek 
to determine planet topography given two images taken from two different locations 
with two different lighting conditions. The algorithms each employ a single cost 
function that combines the computer vision methods of shape-from-shading and stereo 
in different ways. The algorithms are closely coupled and take into account all the 
constraints of the photo-topography problem. 

One such algorithm, the z-only algorithm, can accurately and robustly estimate 
the height of a surface from two given images. Results of running the algorithms on 
four synthetic test image sets of varying difficulty are presented. 
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We know of no universe, whatsoever, 

where music is not in the center of the rocks. 

— Brian Swimme, Canticle to the Cosmos (1990) 
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Chapter 1 
Introduction 



Over the years computer vision researchers have developed algorithms to interpret 
gray scale and color images. These algorithms can be used to extract various types 
of information from the images such as relative motion between images (motion- 
vision), depth from a pair of images (binocular stereo), surface orientation from a 
set of images with different lighting (photometric stereo), surface orientation from a 
single image (shape-from-shading and other shape-from-X algorithms). Most of these 
algorithms work well in specialized circumstances but not as well in general. They 
are not very robust. They work well if everything is perfect, but don't do very well 
in sub-optimal conditions. They operate with many restrictive assumptions, which 
limit their applicability to many real problems. A way is needed to create more robust 
methods that can be applied to wide range of input images. 

One obvious approach is to research ways of enhancing the existing methods to 
make them applicable to a wider range of inputs. Another approach is to combine one 
or more methods to produce a hybrid method with better, more robust performance. 
This thesis looks into the latter approach. The resulting algorithms are called fused 
algorithms. Fusing of two or more computer vision methods, if they bring different 
information to the problem, can create a more robust solution. It's like Kalman 
filtering: with more sensors you get a better estimate. 

This thesis explores and answers the following questions: 

1. What are the different fusion paradigms and how do they differ? 

2. Which fusion paradigm are best for vision problems? 

3. How much additional performance and robustness can be obtained via fusion? 

There are two main approaches to fusing disparate algorithms. I call them the 
module-based approach and the variational approach. In the module-based approach, 
the fused algorithms are run separately on the input images. The various outputs, 
(they might be a sparse depth map from binocular stereo and a dense surface ori- 
entation from shape-from-shading), are then combined somehow to generate a single 
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Figure 1-1: Generic module-based fusion method flow chart. 
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Figure 1-2: Generic variational-based fusion method flow chart. 



solution (see Figure 1-1). The module-based approach is easy to implement since 
existing algorithms can be put together in an ad-hoc manner without having to hnow 
much about the inner workings of each algorithm. The outputs are then combined 
using physically derived constraints to generate the desired combined solution. The 
approach can be applied successfully, but has the disadvantage that this approach 
doesn't fully exploit the information coupling between the methods. 

The variational approach, on the other hand, closely couples the algorithms to- 
gether (see Figure 1-2). This is achieved by formulating a combined cost function 
based on the cost functions of the separate algorithms to be fused. The result is a 
combined optimization problem which takes into account both the explicit and im- 
plicit constraints between the methods. Variational methods, by their nature, can 
exploit any orthogonalities (information content) in the methods. By exploiting all 
the information, the variational approach has the potential to create robust, well per- 
forming combinations of algorithms which can be applied to a wide range of input 
images. How to create such algorithms is the focus of this thesis. 



1.1. PHOTO-TOPOGRAPHY 
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Geometry for left image. 



Geometry for right image. 



Figure 1-3: Example camera and light source geometry for the two images of photo- 
topography. 



1.1 Photo-topography 

A problem that can benefit greatly from a fused vision algorithm is the problem of 
photo-topography. Photo-topography seehs to determine the topography of a planet's 
surface based on two images of the planet, tahen from two different vantage points at 
two different times. This situation is illustrated in Figure 1-3. 

In some ways, this situation is analogous to binocular stereo. Unfortunately, the 
images are typically tahen with two different light source positions as shown in the 
figure. The two images that result will, in general, looh quite different from each 
other even though they are images of the same surface patch. One such set of images 
is shown in Figure 9-1 on page 121. 

Contrast this with the situation that is normally true for binocular stereo images. 
Stereo image pairs are usually tahen simultaneously (or nearly simultaneously), from 
positions that are near to each other, and with the same lighting. The images that 
result looh very similar to each other except for the relative shift of objects (i.e, the 
disparity) due to their distance from the cameras. If the disparity for all points in 
the images and the relative geometry of the cameras is hnown, then the depth of the 
objects can be computed directly. The images are similar, so most stereo algorithms 
determine the disparity by trying to match features in one image with features in the 
other. Since the images from photo-topography looh different from each other, this 
matching approach doesn't worh. 

Photo-topography also shares some aspects with a shape-from- shading problem. 
Shape- from-shading tahes a gray-scale image of a surface and determines the surface 
topography by exploiting the shading information in the image. Shape-from-shading 
requires that the surface reflectance properties be hnown. Assuming the reflectance 
properties of the planet's surface are hnown, we could use a shape-from-shading algo- 
rithm to estimate the surface topography from each of the photo-topography images. 
Unfortunately, the surface estimates based on each image will typically be different. 
They may not even be very similar. In the worse case, the surface estimates from 
each image may not have the same orientation; one could be concave while the other 
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is convex. Thus planet photo-topography is a perfect problem on which to test a 
fusion algorithm based on shape-from-shading and stereo. 

Photo-topography also has aspects in common with a photometric stereo prob- 
lem. Photometric stereo uses two images of the same scene taken with two different 
lighting conditions. Both images are from the same vantage point; the camera is not 
moved between images. The result is two images that look different but where the 
correspondence is explicitly known (a position in one image is always matched with 
the same position in the other image). If the light positions are far enough apart, it 
is possible to determine the surface orientation directly. For Lambertian reflectance, 
two images can constrain the surface orientation to two possible values at each point. 
With three images, it is possible to find a unique surface orientation at each point. 

The photo-topography images have two different light source positions, but the 
correspondence is based on binocular stereo. Like photometric stereo, the two light 
sources can constrain the surface orientation, but only if the correspondence is known. 
Thus in photo-topography, photometric stereo and binocular stereo are closely linked. 

Currently, photo-topography problems are solved using a characteristic strips 
method [Davis and Soderblom, 1984] or more modern shape-from-shading meth- 
ods [Van Hove and Carlotto, 1986]. The characteristic strips method uses a general 
method of solving partial differential equations and is similar to some early shape- 
from-shading algorithms. In this method, the solution is found along strips starting 
from a given (or known) point. The solutions obtained are local, and may not cover 
the space of interest. 

1.2 Related Research 

This research is most closely related to the work of Horn [ikeuchi and Horn, 1981], 
[Horn and Brooks, 1986], [Horn, 1986], [Horn, 1989], Gennert [Gennert, 1987], and 
Szeliski [Szeliski, 1990], [Szeliski, 1991]. The variational (least squares) approach I 
have taken is based significantly on the work of Horn [Horn, 1989], [Horn and Schunk, 
1981], [Negahdaripour and Horn, 1985] and on insights gained from my background in 
control/estimation. The shape-from-shading part of this thesis builds upon the work 
of Horn [Horn, 1989], Szeliski [Szeliski, 1991], and Leclerc and Bobick [Leclerc and 
Bobick, 1991]. The stereo part of this thesis builds on the gray-scale stereo algorithm 
of Gennert [Gennert, 1987]. I use the hierarchical basis functions of Szeliski and use 
conjugate gradient optimization, as do Leclerc and Bobick [Leclerc and Bobick, 1991]. 

This research is also related to the work of Hartt and Carlotto [Hartt and Carlotto, 
1989], [Heipke, 1992], Wildey [Wildey, 1973], and McEwen [McEwen, 1985] in that 
they try to solve the same problem; that is, the determination of planet topography 
from a pair of images. 

Hartt and Carlotto use a Bayesian formulation with multiple images that in some 
ways is very similar to the methods I use. Their cost function and their approach is 
however very different. The Bayesian formulation leads them to a cost function of 



1.2. RELATED RESEARCH 17 

the form (using my notation) 

U ( z ) = II £ \E {k \x, y) - Rk(p, q)\ + \\Vz(x, y)\ 2 dx dy. (1.1) 

JJ k 

While they discuss methods of estimating the albedo of the surface, results are shown 
for constant albedo surfaces with Lambertian reflectance. 1 

They solve their problem using a Markov Random Field Model via a Gibbs distri- 
bution [Geman and Geman, 1984]. They use the Metropolis optimization algorithm 
[Metropolis et al., 1953] with a coarse-to-fine multi-resolution strategy to obtain their 
results. They show results for the algorithm on multiple (2) images based on a digital 
elevation map of upstate New York. They obtain better estimates when two images 
are used instead of one. It is hard to determine from the results they show, but it 
appears that their estimates contain quite a bit of error (only the estimated images 
are shown). The efficiency of their algorithm is also hard to gauge. They state in 
their paper that approximately. 200 sweeps of each image is performed at each reso- 
lution. Inferring that they used 60-by-60 pixel images, with 3 resolution levels, and 
each sweep evaluated the cost function n times where n is the number of pixels on 
that level, the number of function evaluations they required is approximately 

iV^200*(60 2 + 30 2 + 15 2 ) (1.2) 

ss 945000 function evaluations. (1-3) 

This is significantly more function evaluations that the algorithms I have developed 
(my algorithms typically require at most 1500 function evaluations for 65-by-65 pixel 
images and much less on some images). 

Heipke discusses the solution to a multiple image photo-topography problem. In 
the paper [Heipke, 1992], he states that he uses a least squares cost function. In 
addition, the constraint equations that he develops are the same as mine except that 
they apply to more than two images. Since Heipke is not very explicit in this paper, it 
is difficult to tell exactly how he solves the problem. His results are encouraging but 
also show that his method is very sensitive to noise in the images. Heipke's algorithm 
shares many things in common with the algorithms I develop in this thesis since he 
starts with the same constraint equations. Even so, his algorithm doesn't perform as 
well as my algorithms. 

The work of McEwen is based on traditional photoclinometry methods and is 
thus profile based. In [McEwen, 1985], McEwen uses information from two images of 
the same location under different lighting conditions to distinguish between albedo 
and reflectance effects. The image sets McEwen uses are taken from nearly the same 
location so that the pixel correspondence is easy. In many ways, McEwen's methods 
parallel photometric stereo much more than my algorithms. 



1 I make similar assumptions for most of this thesis. 
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Wildey basically solves a stereo problem [Wildey, 1973]. He assumes that that 
the light source is in the same position for each image and tries to determine pixel 
correspondence by matching image brightness. His method estimates the topography 
of the planet surface as it proceeds. He mentions an enhancement to this method 
that would take into account the estimated brightness of the reconstructed surface. 
This enhancement is closely related to my algorithms (and shape-from-shading) and 
would be applicable to images with similar lighting. 

Another related paper is the shape-from-shading and stereo fusion algorithm of 
Grimson [Grimson, 1984]. In that paper, Grimson describes a method for enhancing 
the surface reconstruction step required after feature-based stereo using shading in- 
formation. His method uses the shading information only along features to constrain 
the surface reconstruction and is loosely based on the methods of photometric stereo 
[Woodham, 1980]. 

1.2.1 Relation of this thesis to sensor fusion 

The goal of this thesis is to investigate ways of combining or fusing two different 
computer vision methods. Thus this thesis is related to sensor fusion techniques (also 
called data fusion or information fusion in some contexts). Sensor fusion is a way of 
combining the information from multiple sensors to create better parameter estimates 
than can be achieved with the individual sensors alone. 

Fusion has a much longer history than its use in vision algorithms. Researchers 
and engineers have been combining disparate sensors to obtain better estimates of 
the outside world for a long time. In the literature, three different approaches can 
be identified; 1) Fusion via estimation theory, 2) Fusion via decision theory, and 2) 
fusion via artificial intelligence methods. 

The traditional approach involves the use of Kalman Filtering [Gelb, 1974]. In 
Kalman filtering, or estimation as it more broadly called, the sensor physics are an- 
alyzed and a noise model is postulated. Based on this information, a cost function 
is formulated. Typical cost functions are based on maximum likelihood estimation, 
maximum a priori estimation, minimum variance estimation, or least squares estima- 
tion. For simple noise models (such as Gaussian white noise), these cost functions 
have the form, 

minj= f L(x,...) (1.4) 



possibly subject to additional constraints imposed by the physics. L(x } . . .) is some 
non-linear functional. This same type of cost function is formed when solving vision 
problems using the variational approach [Horn, 1989]. Estimation theory-based fusion 
includes methods that rely on Bayesian statistics [Lee, 1990], [Richardson and Marsh, 
1988], [Hung et al., 1988], as well as least-squares approaches [Shaw et al., 1988]. 

Decision theory approaches use Bayesian Reasoning [Thomopoulos, 1989], Shafer- 
Dempster Reasoning [Bogler, 1987] or ad hoc methods [Mitiche and Aggarwal, 1986]. 
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The artificial intelligence methods range from rule-based hypotheses [Belknap et al., 
1986], [Belknap et al., 1986], [Nandhakumar and Aggarwal, 1988], to blackboard 
schemes [Harmon et al, 1986], to knowledge representation approaches [Pau, 1989], 
to artificial neural networks [Brown et al., 1991]. Sensor fusion is an important 
problem for mobile robots and military threat detection as evidenced by the above 
references. 

An important problem in sensor fusion that is not addressed in this thesis is the 
fusion of widely differing sensors. The problem is particularly difficult for mobile 
robots where the sensors may provide information at several conceptual levels. For 
example, sensors may be available for wheel rotation, range images, visual images, 
touch and sound. Combining information from such disparate sources is much more 
difficult than combining information from sources that operate on the same level (such 
as combining infrared, range and visual images). 

1.2.2 Relation of this thesis to vision fusion schemes 

The computer vision literature is roughly segmented along module boundaries. These 
modules are related to our current ideas about how human vision works. For in- 
stance there are algorithms for binocular stereo, motion vision, and a whole host 
of shape-from-X algorithms. Each of these algorithms performs reasonably well in 
ideal conditions but can degrade quickly in typical real-world situations. It has been 
suggested by many researchers [Aloimonos and Basu, 1988], [Waxman and Duncan, 
1986] that a combination of these modules will do a better job of estimating the ex- 
ternal world. Of course, these ideas are firmly supported by estimation theory where 
it can be easily shown that adding sensors cannot degrade an estimate and most often 
enhances it. 

Along these lines, some researchers in the computer vision community have been 
looking into ways of combining more than one vision cue (and associated algorithm) 
in order to obtain either better or more robust estimates of the external world. The 
techniques for combining the methods range from those that use Bayesian estima- 
tion theory [Matthies and Elfes, 1988], [Hartt and Carlotto, 1989], to those that use 
module-based methods [Moerdler and Kender, 1987], [Moerdler and Boult, 1988], 
[Crimson, 1984]. Another common technique is to use an analytical approach [Wax- 
man and Duncan, 1986], [Aloimonos and Basu, 1988], [Hu and Shrikhande, 1990], 
relying on the constraints from the fused cues (under particular assumptions) to gen- 
erate either a unique or a finite number of possible solutions. 

I prefer the estimation-based methods, since all the information available can be 
exploited, and the assumptions behind a particular method can be quantified. These 
methods also rest on a hrm foundation. The module-based methods, in contrast, 
are more ad hoc. The analytical approach can be used within the estimation-based 
approach to constrain the solution. 
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1.2.3 Relation of this thesis to photoclinometry 

One of the reasons for this research is to improve methods of planetary mapping. In 
particular, this research aims to widen the applicability of computer methods to this 
difficult task. Current techniques for planetary mapping work either on symmetric 
objects or along ridges [Davis and Soderblom, 1984], or use the same methods as 
shape- from-shading [Wildey, 1975], [Van Hove and Carlotto, 1986], [Kirk, 1987] al- 
though the theory was developed independently within the geophysics community. 
In that community, shape-from- shading is called two-dimensional photoclinometry. 
Important aspects for photoclinometry are the determination of the true reflectance 
function for a planet [Davis and McEwen, 1984], [Wilson et al., 1985], [McEwen, 
1991], [Helfenstein et al., 1991] and how to deal with non-constant albedo [Davis and 
Soderblom, 1984], [Helfenstein et al, 1991]. 

While there is no doubt that using the correct reflectance function and accounting 
for non-constant albedo is important, most of the algorithms in this thesis are for the 
simplified case of Lambertian reflectance and constant albedo. Both of these simpli- 
fications are not fatal: the algorithms can be easily generalized to other reflectance 
functions and the non-constant albedo case is dealt with in Section 8.1. In fact, the 
algorithms developed will work, without change, with any reflectance function that 
is smooth and that doesn't contain multiple extrema. 



Chapter 2 
Background 



2.1 Coordinate Systems 

The solution to a photo-topography problem is a description of the surface topog- 
raphy. The most straightforward description for the surface is to represent it as a 
height function over some suitable 2-D domain, 1 that is 



z{x,y) 



(2.V 



Within the vision literature, there are two different choices for the 2-D domain. The 
image- centered domain uses the coordinates of the image as the fundamental do- 
main and assigns a depth (or height) value to the surface point that projects to each 
image position. If orthographic projection is used, then the projection mapping is 
straightforward. However if perspective projection is used, then this mapping can be- 
come quite complex (especially for surface normal calculations). The object-centered 
domain uses a coordinate system associated with the object and assigns a surface 



^^Also known as a Monge patch. 
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Figure 2-1: Coordinate system choices. 
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height value to each point in the domain (see Figure 2-1). Projecting points on the 
object to points in the image is straightforward using this representation for either 
orthographic or perspective projection. However, the projected points, won't in gen- 
eral map to the center of each pixel. To obtain values at each pixel then requires 
some type of interpolation. 

In a stereo system, there is also the choice of whether the coordinate system should 
favor one image or the other. Since there is no apriori reason to believe that one image 
has better information than the other image, it is best to choose a neutral coordinate 
system rather than risk biasing the result toward one or the other image. 

2.2 Image Generation Process 

The photo-topography problem is basically an inverse problem. We seek to determine 
the topography of the object that created the images at hand. Those images are based 
on the object (a forward problem). That is, the interaction of light with an object, 
as seen by the viewer creates the image. The physics behind this process are detailed 
in the sections that follow. 

The image generation process can be neatly broken into four stages (see also 
Figure 2-2), 

1. light falls upon an object (object irradiance) 

2. the light interacts with the object and is re-emitted or reflected (object radiance) 

3. the light then travels to the viewer where it is projected onto the image plane, 
(image projection and image irradiance) 

4. the light is absorbed by the material of the image plane and converted into some 
signal that can be sensed (image transduction). 



2.2.1 Object irradiance 

The light that falls on a particular patch of an object depends on the properties of the 
light sources that are visible (i.e., unobstructed) from that patch. Among possible 
light sources are point sources, distributed sources, and other surface patches (e.g., 
interflection). At each point, £, on the surface, and for each direction, s, the irradiance 
distribution function, -E(£; s) [radiance/solid angle] captures this information. For 
instance, the distribution function for a single point source, ignoring interflection, is 
given by 

E(C;s) = X6(s } s ) (2.2) 
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Figure 2-2: Image generation process steps. 

where A is the light intensity, s is the unit vector in the light source direction, and 
8 is a direction vector version of the Dirac delta 2 



6(t, t ) = 0, when t ^ t , 
/(t)$(t,t )rft = /(t ) 



(2.31 



where the integral is over the surface of the sphere S. 



2.2.2 Object Radiance 

The light that is emitted by a surface depends both on the light that impinges on 
the surface and the surface reflectance properties. The irradiance distribution func- 
tion, -E(£; s) captures all the needed information about the light source, while the 
surface reflectance properties can be described using the Bidirectional Reflectance 
Distribution Function (BRDF). The BRDF is an intrinsic function of the surface and 
doesn't depend on surface irradiance. At the point £, the BRDF /(£; n; v; s) relates 
the brightness of the surface patch with normal n illuminated from the direction s 
and as seen from the direction v. Using these two distribution functions, the surface 
radiance [power/solid angle/area] in the direction v can then be defined as 



L(i- n; v) = II f(£; n; v; s)E(£; s)(s • n)^(s) 

JJH(h) 



(2.4) 



2 This delta function is defined on the surface of the sphere. One way to define it is to repre- 
sent the direction t in spherical coordinates, t = (sin 6 cos <f>, sin #sin <f>, cos 0) T , so that we require 
J-7r/2/o' r HH'f' ^) ^o((f>o, o)) cos 9 dcf)d9 = 1, hence, 



6(t(<f>, 6)M<f>oM) = ^ ~ 4>o)6(0 - o )/cos(0 o ), 
where 8(.) is the normal scalar Dirac delta. 
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where H(n) is the hemisphere of possible light source directions for the patch at 
surface point £ with normal n, and dio(s) is the solid angle subtended in the direction 
s. Representing s in spherical coordinates this expression becomes 

L((;h;v) = T T /(£; n; v; s)£(£; s) sin # cos 6d6d<t>. (2.5) 

J--K JO 

The two most common reflectance models are the Lambertian model for matte 
surfaces, and the specular model for shiny metallic surfaces. In the Lambertian case, 
the BRDF is independent of the light source direction, normal direction and viewing 
direction, 

f(C;h;v;s) = -p(0 (2.6) 

7T 

where /?(£) is the surface albedo, or the fraction of light re-emitted by the surface. 
Evaluating Equation 2.5 for a Lambertian surface illuminated by a single point source 
at infinity, the surface radiance is found to be 

I((;n;v) = ^(s .n). (2.7) 

7T 

In the specular case, all the light from the direction s is reflected into the direction 
2(n • s )n — s so the BRDF is 

f(£; n; v; s) = />(£)£( v, 2(n • s )n - s ). (2.8) 

For a single point source at inhnity, the surface radiance for a such a surface is found 
to be 

L((; n; v) = A/>(£)*(v, 2(n • s )n - s ) (2.9) 

In general, the albedo of a Lambertian surface is different from the albedo of a specular 
surface. 

See [Wilson et a/., 1985], [McEwen, 1985], and [Davis and McEwen, 1984] for a 
sampling of the types of radiance functions used within the photoclinometry held. 

2.2.3 Reflectance Map 

The preceding representation of the surface radiance is a local representation in that 
the radiance is defined based on local surface orientation, viewer direction and light 
source direction. Many vision researchers have found it convenient to use a represen- 
tation of the surface radiance based on a global coordinate system. This alternate 
representation is called the Reflectance function. Given known surface properties and 
a known light source distribution, the reflectance function, R(£, n; v), is defined using 
the corresponding surface radiance via a change in coordinates, 

R(t G ;n G ;v G ) = mn;v) (2.10) 
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Figure 2-3: Two Perspective projection geometries. 

where £q, iig, and Vq are the surface position, normal vector, and viewing direction 
in global coordinates. 

2.2.4 Image Projection 

Points in the image plane are related to points on the object via perspective projection. 
Figure 2-3 shows two perspective projection geometries. Figure 2-3(a) shows the 
projection geometry for a camera where the image plane is behind the lens. The 
projection of the center of the lens into the image plane is called the principal point. 3 
Choosing the origin of the image coordinate system to be at the principal point leads 
to the simple equations presented in this section. 4 Note that objects are projected 
onto the image plane in an upside down orientation. 

An equivalent geometry is shown in Figure 2-3(b). The only difference between 
the camera geometry shown in Figure 2-3(a) and this geometry is that the image plane 
has been moved to in front of the lens. While this geometry is not physically possible, 



3 The center of projection is at the back nodal point for thick lenses. 

4 Determining the position of the principal point for a given camera is part of the classic interior 



orientation problem (see [Horn, 1986]). 
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it produces exactly the same image (as can be seen by similar triangles) except that 
the image is no longer upside down. Mathematically, the difference between these two 
geometries is just a sign change in the equations. For instance, the image position r 
of the object point R in Figure 2-3(a) is 

r R 

— = — , (2-11) 

r-z R z' v ; 

while it is 

r R 

— = — (2.12) 

r-z R-z v ; 

for the geometry in Figure 2-3(b). These are Perspective Projection Equations. In 
both expressions, z is a unit vector in the direction of the negative optical axis. In 
the equations that follow, the geometry of Figure 2-3(b) will be used to avoid having 
to keep track of the minus sign in the perspective projection equation. 

The equations above can be simplified if a special coordinate system is used, 
namely a coordinate system with origin at the principal point, and orientation such 
that the x and y axes are aligned with the image coordinates and with the negative 
z axis along the optical axis of the lens. This coordinate system is evident in the 
labels assigned to the axes in Figure 2-3 (b). For this coordinate system, r • z = / (a 
negative number). Expanded, these equations are 

x/f \ (X/Z\ 

y/f \ = \ Y / z (2-i3) 

/// / V z/z ) 

when r = (x, y, f) T and R = (X, Y, Z) T . 
2.2.5 Orthographic Projection 

The projection equations can be further simplified when the depth range of the object 
is small compared to the distance of the object from the camera. The resulting 
approximation is called Orthographic Projection. Consider the hrst order Taylor's 
series expansion of the perspective projection equation about a nominal depth Z 0} 

7 = ^ + ^(Z - Z ) + 0((Z - Z ) 2 ) (2.14) 

where R = (X, Y, Z ) T . When (Z — Z ) <C Z , the hrst order term can be neglected 
resulting in the orthographic approximation, 

r R 

7 = Y <™> 
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Figure 2-4: Image irradiance. 

Orthographic projection is convenient for computer vision problems since the map- 
ping is linear. In fact, most shape- from-shading algorithms assume orthographic 
projection. Unfortunately, the full perspective projection must be used for photo- 
topography to avoid throwing out the very information we wish to estimate (i.e., the 
stereo depth information). 

2.2.6 Image Irradiance 

The mapping between object points and image points is only half of the story for 
image generation. We also need to know how the brightness of the image is affected 
by the lens. Figure 2-4 shows that an object patch of area 80 projects to an image 
patch of area 81. Assuming a perfect lens, the image irradiance [power/area] from 
this patch is 



£(r)=i^(r);n;v(r)) 



7T 



d 



H/, 



cos 4 a 



(2.161 



where _L(£(r); n; v(r)) is the radiance of the corresponding object patch, d is the 
diameter of the lens, and a is the off-axis angle of the projecting ray. In this equation, 
v and £ are functions of r via the projection from image points r to object points. 
Equivalently, this expression can be restated using the reflectance function rather 
than the radiance function, 



£(r) = i?(C(r);n;v(r)) 



it J d 

H7, 



cos 4 a 



(2.17) 



2.2.7 Image Transduction 

The final stage in producing an image is the conversion of light into a signal that can 
be used, namely digital form. Whether the digital images are scanned photographs 
or are obtained directly via a digital camera, they will contain distortion (either 
spatially or in color). The effect of this distortion can be removed via calibration 
so that the measured image irradiance can be related to the object radiance in a 
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Figure 2-5: Stereo geometry. 

straightforward way. With perfect calibration, the image irradiance equation can be 
put in the especially simple form, 



£(r) = i?(C(r);n;v(r)) 



(2.181 



where the cos 4 a and constant scale factors have been removed as part of the calibra- 
tion. Equation 2.18 is referred to as the Image Irradiance Equation. 

2.3 Stereo 

The images in a photo-topography problem are taken from two different vantage 
points (see Figure 2-5). If the relative position of the cameras is known, and it is 
known which pixels in the left image correspond to which pixels in the right image, 
then it is possible to determine the depth directly for the surface points that map 
to those pixels. This is the basis of binocular stereo. Of course, the hard part is 
determining the correspondence between pixels. 

For normal stereo situations, the cameras are close together and both pictures 
are taken simultaneously. 5 The stereo images that result look very similar, mostly 
differing in the shift of objects in each image caused by perspective projection. The 
difference in the shift of an object point in the left image and the right image is called 
the disparity. The photo-topography images are taken with cameras that are often far 



3 The stereo images of aerial photography are taken nearly simultaneously. 



2.3. STEREO 29 

apart and at different times. In this case, it is possible (because of differing lighting) 
for the two images to look very different. This makes the correspondence problem 
even harder. 

To see how the depth can be determined directly from the disparity consider 
Figure 2-5. Points in each image map to rays in 3-space. If corresponding points 
in each image map to rays that intersect, the depth can be determined using simple 
geometry. If the rays don't intersect, we can use geometry to determine the depth of 
the point that is halfway between the rays at their closest approach. 

First, find the relationship between the two camera coordinate systems. Suppose 
we know the position of the principal point of each camera in some global coordinate 
system, Pi and P 2 , and we also know the rotational transformation matrices from 
each local camera coordinate system to the global coordinate system, T\ and T 2 , then 
the coordinates of the point, £, in the two camera coordinate systems is 

R^rr^-Pi) 

R 2 = r- 1 (e-P 2 ). (2 - 19) 

Equations 2.19 are the Stereo Constraint Equations which relate the position of points 
in the global coordinate system to points in each of the local camera coordinate 
systems. 

By removing £ from the above equations and defining b = P 2 — P 1? the relationship 
between a point in Camera Coordinate System I and a point in Camera Coordinate 
System 2 is found to be 

R x = T^b + T 1 - 1 r 2 R 2 . (2.20) 

Now, determine the relationship between disparity and depth. Suppose we are 
given image points, r 1 = (x 1 ,y 1 , f) T and r 2 = (s 2 ,J/ 2 ,/) T (one in each image), that 
correspond to the same surface point, then the best estimate for the surface position 
can be found by finding the point on each ray (along r x and r 2 ) where the distance 
between the rays is minimized. That is, the problem 

minll-sri + b + tr 2 || 2 (2.21) 

s,t 

must be solved where s and t are scalar parameters. For now assume all the vectors 
are given on the same basis. 

By differentiating the above equation with respect to s and t, setting the resulting 
equations to zero, and solving, it is found that the minimum occurs when 



(r 2 • r 2 )(b • ri) - (r x • r 2 )(b • r 2 ) (r 2 x b) • (r 2 x r x ) 



(ri-ri)(r 2 -r 2 ) -(ri -r 2 ) 2 ||r 2 x r x || 2 

(r x • r 2 )(b • ri) - (r x • r x )(b • r 2 ) (r x x b) • (r 2 x r x ) 



(2.221 



(r x • ri)(r 2 • r 2 ) - (r x • r 2 ) 2 ||r 2 x r x 
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If i - !, r 2 , and b are coplanar then the above set of equations is just a fancy way of 
writing the Law of Sines. 

The global position of the point halfway between the two rays at the closest 
approach is then 

£«P 1 + sr 1 + -(-sr 1 + b + *r 2 ), (2.23) 



b 1 (T 2 r 2 x b) • (T 2 r 2 x T iri ) 
2 2 ||T 2 r 2 x T iri | 

i (r iri x b) • (r 2 r 2 x r iri ; 

2 ||T 2 r2 x TirJ 



Pi + " + ^ i i :s j v i i ; ' i; r iri 



+ o ^ ™ — wi ^2r 2 (2.24) 



where Equation 2.24 has been written for P 1? P 2 , and b given in global coordinates, 
and i - !, r 2 given in the appropriate local camera coordinate system. 

The equations simplify greatly if the r 1? r 2 , and b are coplanar (i.e., if (r 2 X r x ) b = 
0), the cameras are aligned so that the optical axes are in the direction — z, and x is 
along b (y is chosen to complete the right-handed coordinate system). This is exactly 
the stereo geometry that is assumed to exist for most binocular stereo algorithms. 
With these restrictions, the equation above becomes 

f = Pl+ b + |£l±£4. (2,25) 

2 2(x 1 - x 2 ) 

where r 1 = (x 1 , y 1 , f) T , r 2 = (x 2 , j/ 2 , /) T , and b = (6, 0,0) T . The quantity (x 2 — Xi) 
in the above equation is the disparity mentioned earlier. In this form, the depth z is 
found to be 

, = £ . i = Pl . S+ » + _^_ y . (,,26) 

Note that the disparity can be mapped directly into depth only in this special situa- 
tion. For more general situations, Equation 2.24 must be used. 

2.4 Photo-topography Problem Formulaton 

We now have enough background to formulate the photo-topography problem. The 
problem to be solved is: 

Given two images of an area on the planet's surface, taken at two differ- 
ent times from two different positions, determine the topography of that 
section of the planet's surface. 

The solution is constrained by geometry and the image generation process. Specifi- 
cally each image is constrained by the perspective projection equation (Equation 2.12), 
the image irradiance equation (Equation 2.18) and the stereo constraint equations 
(Equations 2.19). 
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Combining these equations we find that the photo-topography problem is constr- 
ained such that 

E^(r 1 ) = R^(t;n;-T 1 r 1 ) 

E^(r 2 ) = R^(C;n;-T 2 r 2 ) (2 " 27) 



where 



Tl 



r 2 



f Tl -\C-Pi)_ 

* i 

2.28 
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2 )-r 1 z 2 ' 



are the local image plane position vectors (one in each camera coordinate system), 
E^ 1 ' and E^ 2 > are the image brightness measured in the hrst and second cameras 
respectively, and R^ 1 ' and R^ 2 ' are the reflectance maps based on the hrst and second 
light source positions. As before, £, is the surface position in global coordinates. 

Equations 2.27 and 2.28 are the General Photo-topography Constraint equations. 
They state that in the absence of noise, the photo-topography images E^ 1 ' are created 
by the interaction of light with the underlying surface z(£) via the reflectance maps 
RS 1 ' and as seen by the two cameras. Note that these equations can be used to 
generate a set of photo-topography images given the surface description, z(£). In 
contrast, determining the surface topography given the images, camera geometry and 
reflectance properties, is an inverse problem. 

It's instructive to review the assumptions behind these equations. The perspec- 
tive projection equations assume perfect lenses and perfect knowledge of the camera 
principal points and optical axes. The surface radiance equation assumes we have 
perfect knowledge of the surface reflectance properties, light source directions, and 
all surface points are visible from both cameras (i.e., there are no self occlusions). 
The simplified form of the image irradiance equation assumes we either have a per- 
fect sensor or we can perfectly calibrate the sensor to remove any abnormalities from 
the sensor/lens combination. The stereo equations assume we know the relative posi- 
tion and orientation of the two cameras perfectly. The only assumption that is truly 
artificial is the assumption of perfect knowledge of the reflectance maps. 6 With more 
careful measurements and more expensive equipment it is possible to approach per- 
fect knowledge of the other assumptions. The assumption that all surface points be 
visible merely restricts the bumpiness of the surface that this research is applicable 
to. 



6 Especially since the effects of non-uniform albedo and interflection are buried in the reflectance 
maps. 
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Figure 2-6: Introduced global coordinate system. 

2.5 Simplifications 

The equations in the last section dealt with the general case (subject to the assump- 
tions mentioned in that section). There are several simplifications that make the 
equations easier to solve. 

2.5.1 Special global coordinate system. 

So far all the equations have been written for any global coordinate system. I would 
now like to restrict the equations to a particular global coordinate system, namely the 
coordinate system shown in Figure 2-6. This coordinate system is defined as follows 

I. Place the origin of the global coordinate system half way between the principal 
points of the two cameras. (That is, place the origin at the point P x + b/2.) 



2. Choose the x direction along the line connecting the two cameras, 



x 



b/||b| 



(2.291 



3. Choose z as the average optical axis direction of the two cameras projected 
into the plane perpendicular to x . (That is, 



z x • x )x ) + (z 2 



z 2 • x )x 



Z X • XoJXo) + (z 2 



(2.301 



z 2 • x )x 



4. Choose y in the direction of z X x in order to create a right handed coordinate 
system. 

5. Also set up a virtual image plane with / = 1. 

In this coordinate system, Pi = —b/2, P 2 = b/2, and b = (6, 0,0) T . 
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2.5.2 Removing the view direction dependence. 

A common simplification for computer vision is the assumption of a Lambertian 
reflectance map. Since a Lambertian surface reflects light equally in all directions, 
we see from Equation 2.7 that the radiance function is not dependent on the viewing 
direction. Thus the dependence on v can be removed from all the equations. 7 

Note that the viewing direction can also be removed from the equations when the 
held of view is small. In this case, the viewing direction is approximately the same 
for all points on the surface and its effect can be subsumed into the reflectance map. 
Doing this would, of course, introduce an error into the calculations. This error would 
be small for photo-topography since the cameras are so far away from the surface. 
The large viewing distance requires the use of a telephoto lens which has a small field 
of view. 

2.5.3 Constant albedo. 

Thus far the equations have included terms that denote position on the surface £. 
The main reason for this dependence is to take into account varying albedo, varying 
reflectance properties, or both. To simplify the situation we could assume that the 
reflectance properties, albedo, or both are constant across the surface. Assuming the 
reflectance properties, but not the albedo, are constant across the surface results in 
a reflectance function that is separated, 

R(C;n;v) = p(C)R(n;v) (2.31) 

where i?(n; v) is the reflectance function for a surface with uniform albedo, no in- 
terflection, and no self occlusion. As for i?, any light source effects are included in 
R. 

When both the reflectance and albedo are constant, the dependence of the re- 
flectance map on surface position can be removed, 

i?(C;n;v) = i?(n;v). (2.32) 

Combined with either Lambertian reflectance 8 or when the field of view is small, the 
dependence on v can be dropped also, 

R(£;h;v) = R(h). (2.33) 

This is the representation of the reflectance function that is seen most often in the 
vision literature. The simplification restricts the applicability of that research so 

7 This is true for any reflectance function that is view independent, not just Lambertian 
reflectance. 

8 Or any other view independent reflectance function. 
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Figure 2-7: Coordinate systems with aligned optical axes. 

that only uniformly colored surface patches have the possibility of being estimated 
correctly. When algorithms based on this simplification are applied to images that 
violate these simplifications, we would expect errors at the transition between different 
colored parts of the surface, within differently colored areas, or both. Trying to 
estimate both the surface topography and surface albedo substantially increases the 
number of unknown variables (possibly by a factor of 2) and significantly slows down 
the convergence (see Section 8.1). 

2.5.4 Aligned cameras. 

The final simplification that can be made is to align the cameras so that their optical 
axes are parallel (which will also be parallel to the global coordinate system's z axis 
by construction). This coordinate system is shown in Figure 2-7. When the cameras 
are aligned, the rotational transformations T\ and T 2 are identity transformations 
which simplifies the stereo equations (Equations 2.19) considerably to 



Ri=£-Pi, 

R 2 = £ - P 2 



(2.34) 



While this situation is very unrealistic for the photo-topography problem, any set of 
images can be re-projected into this coordinate system (see Section 8.5). Hence this 
simplification does not seriously restrict the applicability of the research. 



2.5.5 The simplified equations. 

The rest of the thesis is based on equations that take into account all of the simplifi- 
cations described. In particular, the simplifications in the following list are made: 
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• A special global coordinate system that is halfway between the two camera 
positions is used. 

• All surface points are assumed to be visible from the two cameras. 

• The reflectance properties of the surface are assumed to be constant with no 
interflection and no mutual occlusion. 

• The radiance properties of the surface are assumed to be Lambertian, allowing 
the view direction dependence to be dropped from the reflectance equations. 

• The surface is assumed to have constant albedo allowing the position depen- 
dence to be dropped from the equations. 

• The camera optical axes are assumed to be aligned with each other allowing the 
rotational transforms to be dropped from the stereo constraint equations. 

Unless noted, the rest of the thesis assumes that all these simplifications hold. I will 
specifically point out results that apply to the more general case. Taking into account 
all of these simplifications we find that the simplified photo-topography problem is 
constrained such that 

£W(n) =i? (1) (n) 

E^(r 2 )=R^(n) (2 - 35) 



where 

r 2 



(£ + b/2)/ C + b/2 



(t-b/2)f C-b/2 [Z -^ 



(Z - b/2) • z 2 i • z 2 
If we define z = £ ■ z and r = /£/z, then the constraint equations can be written 

£ (1 V+^) = i? (1) (A), 

}l (2-37) 

£( 2 V-g) = i?( 2 )(n). 

The definitions for z and r above introduces a peculiarity into the equations: the 
units for z and r are not the same. In particular, z, 6, and £ are in planet units (say 
miles), while r and / are in camera units (say millimeters). 

I find it convenient in subsequent chapters to use a slightly different version of 
these last equations where the components of r = (x, y, f) are explicit and the normal 
vector n is parameterized using gradient components p and q } 

^fp 2 + q 2 + I 



36 



CHAPTER 2. BACKGROUND 



































s :± 





Pixel Coordinates 



Image Plane 1 ' 




Object ; 
Common surface point 



Image Plane 2 



Figure 2-8: Camera calibration geometry. 



where 



P 



fz* 



%z x -\- z 

f z y 
yz y + z' 



f2.39] 



Equation 2.39 are referred to as the Integrability Constraint Equations since they 
relate the surface normal components to the partial derivatives of the surface height 
z. The equations can be used constraint the values of p and q that are consistent 
with an underlying surface. 

Writing Equations 2.37 in terms of (x, y) and p and q produces the Photo-topography 
Constraint Equations 



E {1 \x + 
E {2 Hx- 



2z 
2z 



>y 



>y) 



RfV(p, 

R {2) (p, 



<i) 



? 



(2.40) 



I will be working with these equations in the chapters that follow. 



2.6 Camera Calibration 

In order to relate positions in the image to direction vectors in 3-space, the origin 
of the camera coordinate system must be known. Finding this origin is part of the 
interior orientation problem of classical optics. As mentioned in Section 2.2.4, this 
origin is the principal point in the image plane. In the special coordinate system 
of Figure 2-8, the position of the image plane origin with respect to each camera 
coordinate system origin can be specified by a vector v = (v x , v y , f) T . These vectors 
specify the offset (in pixel coordinates) of the image plane for each camera. Suppose 
u is the pixel position of an object point in the global image and v is the position 
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of the origin of the global image plane in global coordinates, then given the pixel 
position of the projection of this this same object point in the camera images u x and 
u 2 , the offset of each image plane is 

fb 
Vi =v + u -u 1 + -^— , (2.41) 

zZ 

fb 

v 2 = v + u - u 2 - —=-. (2.42) 

The values of v x and v 2 can be quite large in the aligned coordinate system indicating 
that the images must be shifted far away from the camera coordinate system origin. 9 
While this is not possible physically, it is a consequence of re-projecting real images 
into the aligned coordinate system. 10 



9 This is especially true with the field of view is small. 
10 See Section 8.5 for how to do this. 
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Chapter 3 
Fusion Strategy 



This chapter presents my general solution strategy to the photo-topography problem. 
The idea is to closely couple the solution of shape-from-shading and stereo so that 
each can help the other. The important thing is to take into account the strengths of 
each method and use each method to cover the weaknesses of the others. 

The hrst step in solving the photo-topography problem is to analyze the infor- 
mation available as part of the problem. The photo-topography problem has two 
information sources: 1) the gray levels in the each image are an indication of the 
surface orientation with respect to the light source and 2) assuming corresponding 
pixels in each image can be matched up, the stereo information can be used to recover 
the shape. Of the two, the gray level information is more directly accessible since 
no correspondence must be found. Finding the correspondence between pixels is a 
classical problem for stereo algorithms and is very hard. 

The photo-topography images also offer one other source of information. Since 
the images are typically taken with two different light source positions, the gray 
levels of corresponding pixels constrain the set of possible surface orientations as they 
would for photometric stereo [Woodham, 1980]. For Lambertian reflectance, this set 
contains at most two orientations. Since extracting this information relies on pixel 
correspondence, it has the same limitations as the stereo information. 

When designing a fusion algorithm it is also important to analyze whether the 
information sources are independent or not. By independent I mean, 'Is it possible 
to differentiate between the two information sources'? Or said another way, 'Given a 
combined signal containing two or more information sources is it possible to estimate 
the relative contributions of each source'? Or said a third way, 'Given a combined 
signal are the information components observable' ? 

Doing this analysis on the photo-topography images we find that the informa- 
tion sources are indeed independent. The shading information is strongest when the 
shading is smooth, while the stereo information is strongest near surface discontinu- 
ities (which contribute to shading discontinuities) and when the cameras are widely 
separated. The photometric stereo information is strongest when the light source 
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constraint 


Correspondence 
constraint 


Correspondence 
constraint 


Photometric 
stereo 



Table 3-1: Photo-topography problem information sources 

positions are widely separated. 

Given the previous analysis we would expect that the best estimates of the surface 
would be found when the images contain large regions of smooth gray level changes 
intermixed with areas of rapidly changing gray level, the cameras are well separated, 
and the light source positions are separated. 

There are several possible fusion paradigms as I discussed in Chapter 1. The 
most promising fusion paradigms take into account the most amount of information. 
The modularized approaches (where existing vision algorithms are applied to the 
problem individually and where these individual estimates are then combined after- 
wards) exploit some of the information but not all of it. In particular, if you think 
of these problems in terms of a table or matrix as in Table 3-1, the modularized 
approaches exploit the information along the diagonal but do not take into account 
the off-diagonal coupling between methods. On the other hand, I use a close coupled 
approach based on variational calculus that exploits all the coupling inherent in the 
problem. This approach is not a panacea, however, as we will see later. It is still 
possible to choose cost functions that don't encourage cooperation among the var- 
ious information sources. However, by understanding the couplings inherent in the 
problem a suitable cost function can usually be found. 

3.1 Variable representations 

The most important decision when developing a fusion algorithm is the variable 
representation. Choose the right representation and everything will work together 
smoothly. Choose the wrong representation and the algorithm may be hindered by 
slow convergence, local minima, discontinuities, etc. Aside from this performance 
effect, how can a good representation by recognized? I can't offer a general solution, 
but I can offer some rules of thumb. 



• A good representation captures all the necessary information in the problem 
with no redundancy. For example it is better to use a parameterized family of 
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functions rather than a digitized version of the same function if the number of 
parameters would be much less than the number of sampled points. 

• A good representation can easily capture the constraints of the problem. For 
example, it is better to use a function that automatically matches the boundary 
conditions of a problem rather than to impose the constraint using penalty 
functions or Lagrange multipliers. 

• A good representation usually has equations that are simple and easy to under- 
stand. 

• A good representation is as close as possible to the problem. It is not removed 
by several integrations, differentiations, or other mappings that can introduce 
non-linearities, bias, drift, or additional unobservable parameters. 

• A good representation allows for easy information exchange between the differ- 
ent information sources in the problem. Such a representation would directly 
express constraints from one information source so that the other sources will 
automatically take them into account. 

The variable representation defines the space that will be searched. Usually, so- 
lution spaces that are lower dimensional, smoother, and more bowl-like will make 
finding the solution faster and easier. 

I present four different variable representations for the photo-topography problem 
in Chapter 4. The four representations highlight the importance of choosing the right 
representation. 

3.2 Cost functions 

After the choice of variable representation, the choice of cost function is the next 
important. Like choosing the right variable representation, choosing the right cost 
function will affect the performance of the algorithm and the algorithm's robustness. 

Building a closely coupled fusion algorithm, requires that all the governing equa- 
tions, constraints, and desired outcomes be formulated into a single cost function, 
usually some sort of energy-based functional. Since everything is combined together 
to generate a single scalar value (the cost), the algorithm is free to perform trade-offs. 
For instance, the algorithm can trade-off accuracy of the result in order to satisfy 
another constraint (such as smoothness). The types of trade-offs that are allowed 
and/or favored can be controlled by changing the relative weights between terms in 
the cost function, or by totally reformulating the cost function (I show examples of 
both in Chapter 4). 

Cost functions are relatively easy to formulate once a variable representation has 
been chosen. Simply write the governing equations in terms of these variables, and 
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create the cost function directly from these equations. If necessary, sum, integrate or 
somehow combine the results from the constraints to form a scalar cost. It might be 
necessary to add some regularizing terms to the cost function in order to produce a 
well posed problem (see [Tikhonov and Arsenin, 1977] or [Hadamard, 1923]). 

Many vision problems require some type of regularization to constrain the solution. 
Typically, smoothness is assumed for some variable (such as surface height) and is used 
to create a smoothness regularization term. Regularization terms are required when 
there are too many solutions to a given problem. The regularization terms can then 
be used to constrain the set of possible solutions to be physically reasonable. In all the 
photo-topography algorithms presented later in this chapter I have included a surface 
smoothness term. While the smoothness term is not strictly required (since shape- 
from-shading problems do have unique solutions), 1 the regularization term serves to 
speed up convergence by making the solution space more bowl-like. 

Constraints can be added to a candidate cost function using the method of La- 
grange multipliers or by adding a penalty function. The Lagrange multiplier method 
is used for hard constraints, while the penalty functions are used for soft constraints. 
Hard constraints are required to be met exactly by the solution, while soft constraints 
need not be met exactly. Hard constraints restrict the space of the feasible solutions, 
while soft constraints create bowl-like edges on the solution space. Both methods can 
be used to help convergence to a solution and to restrict the search space to more 
desirable solutions. 

Different cost functions can be generated from the same set of equations and 
variable representation by changing the desired outcome. For instance, the cost func- 
tion for determining the surface slopes from the photo-topography problem would be 
different from the cost function for finding the surface height. 

3.3 Solution Techniques 

A typical cost function for a vision problem is of the form 

mmJ=[[L(u;u>;u»;...)dxdy (3.1) 



where u are the optimization variables, and I is a possibly non- linear function of 
the optimization variables and its derivatives. The integral is taken over the domain 
of the optimization variables or the problem. The solution to this problem can be 
found by solving the associated Euler- Lagrange equations (see a variational calculus 
book such as [Courant and Hilbert, 1962] or Horn's Appendix [Horn, 1986] for more 
details). 

The Euler-Lagrange equations for a problem such as the one above are typically 
coupled non-linear equations. Such equations are usually very difficult to solve ana- 

1 When the images contain singular points [Saxberg, 1989]. 
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lytically but can sometimes be solved numerically by converting them into discrete 
equations. The conversion process involves substituting discrete approximations for 
any derivatives of the optimization variables. The optimization variables may have 
to be approximated by a discrete vector as well. The equations are then re-arranged 
to create iterative update equations of the form 

u(k + l) = f(u(k),u(k-l),...) (3.2) 

where u(k) is the value of the optimization variables for the k-th iteration. 

When u has many components and when the components are updated in sequence 
based on the most current estimate u, the resulting update scheme is called a Gauss- 
Seidel optimization. When all of the components of u are updated simultaneously 
based on a previous estimate for u, the resulting scheme is called a Gauss-Jordan 
optimization. Gauss-Seidel optimization schemes have higher convergence rates and 
are more robust, and are best implemented on a serial computer. Gauss- Jordan 
schemes, while they have lower convergence rates and are not as robust, can be 
implemented on parallel computers. Parallel computations can produce results faster. 

Another way of solving the optimization problem posed above is by using direct 
optimization techniques. In this case the cost function, instead of the Euler- Lagrange 
equations, is discretized. Any integrations are approximated by sums and any deriva- 
tives are approximated by differences. The resulting cost function is of the form 



mm 



•7 = ££/(«) (3-3) 



where f(u) is a discrete approximation of L(u; u 1 ; u 1 ; . . .). Any of the wide range of 
optimization algorithms that are in the literature can then be applied to this problem. 
For vision problems, conjugate gradient optimization shows the most promise. The 
conjugate gradient scheme doesn't require the formation of the problem Hessian (a 
linear approximation to the second order derivative of the solution space at a given 
point), which for an optimization problem with N variables is an A-by-A matrix. The 
conjugate gradient scheme is important for vision problems since for a typical 256-by- 
256 image, the shape-from-shading problem would have 256 2 or 65536 optimization 
variables. The Hessian for this problem would have 256 4 or over 4 billion elements! 

Conjugate gradient optimization requires that both the cost function and its gra- 
dient are computable. While an approximation to the gradient of the cost function 
can be computed using finite differences, this approach is usually slow (since it re- 
quires at least A function evaluations for each gradient evaluation). I use analytically 
determined gradients in all the algorithms presented later in this chapter. 2 

One strategy that can also be used to solve non-linear problems is the method of 
homotopy or the continuation method. The continuation method involves solving a 



2 Verified by comparing with numerical gradients. 
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series of related problems of increasing difficulty that have similar solutions. Typically, 
the original problem is restated to include a parameter, such as /3, which controls the 
difficulty of the problem. When (3 = 1 (for instance), the problem is easy to solve 
or has a analytical solution. When (3 = (for instance) the original problem is 
recovered. To solve the original problem, a series of problems are solved while (3 
is slowly decreased. The solution of each subsequent problem is used as the initial 
condition for the next problem. 

One simple way to use the continuation technique is to introduce a regularization 
term based on (3 that creates a convex problem when (3 = 1, 



mm, 



EE/W + ^M- (3-4) 



where g(u) is a regularization function that forces the problem to be convex. 3 During 
optimization, (3 is slowly reduced to zero, where the original problem is recovered. 
The hope with the continuation method is that the solutions to the series of problems 
will be close to the solution of the original problem so that they are good initial 
conditions. 

I use a form of the continuation based method for the algorithms in this thesis. 
Since shape- from-shading problems have a well defined solution when singular points 
are in the image, the smoothness terms are not needed to guarantee a solution. 
The smoothness terms do, however, help convergence. During optimization, I slowly 
reduce the smoothness parameter toward zero in order to avoid biasing the solution 
[Horn and Brooks, 1986] . 

The questions of existence and uniqueness come up when working with optimiza- 
tion algorithms such as those presented in this thesis. For the types of cost functions 
presented in this thesis, it is clear that a solution exists; the cost functions are bounded 
from below by zero. That is, the best possible value for the cost function is zero and 
can be achieved only when the estimated surface images and the actual images match 
exactly and when any regularization terms are set to zero. 

The uniqueness of a solution depends a great deal on the surface to be estimated. 
In general, both global and local minima will exist (as evidenced by the Hard Crater 
problem discussed in Chapter 7). The optimization techniques discussed above only 
guarantee convergence to a local minima. The global minimum may only be achieved 
if the initial conditions for the optimization algorithm are close to the true solution. 

3.4 Speed-up Techniques 

A very well researched part of optimization theory is how to speed up the conver- 
gence. For vision problems there are two promising speed up technologies: the use 

3 Convex problems have the property that the solution space is essentially bowl-like, with no local 
minima and a single global minima which makes them very easy to solve. 
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Figure 3-1: Hierarchical basis functions. Shown are the hierarchical basis nodes at each level 
(circles, triangles, etc.) and the associated interpolation function extent for the front-most 
nodes. The nodal basis for the 9-by-9 domain would have nodes at each pixel location and 
an extent half-again as small as the level 1 extents shown. 



of hierarchical basis functions, and multi-grid methods. Both methods try to speed 
up the optimization problems by increasing the information transfer spatially. The 
methods are based on the property of many vision algorithms that an optimization 
variable within a grid of optimization variables may only be affected by its nearest 
neighbors. Due to the local connectedness, many vision algorithms have diffusion- 
like properties; the solution must diffuse throughout the grid. Schemes that directly 
transfer information over longer distances instead thus may speed up an algorithm. 

Using hierarchical basis functions transforms the optimization space as seen by 
the optimization algorithm but not as seen by the vision algorithm. Basically it's like 
a change of basis. Figure 3-1 shows how a 9-by-9 domain would be represented in hi- 
erarchical basis. In particular, note that the nodes of the hierarchical basis propagate 
information over a much larger range (due to their extended interpolation extent) 
than the nodes in the nodal basis (particularly for the nodes at the upper levels). 
The figure shows linear interpolation between nodes, but any interpolation scheme 
can be used to build a hierarchical basis (see [Szeliski, 1990]). The problem is solved 
as before but with the new variables. If necessary, the variables can be transformed 
to the nodal basis when computing the cost function or gradient. Unfortunately, 
all these transformations have the potential to introduce round-off errors which can 
adversely affect sensitive algorithms. 

The hierarchical basis functions have the most effect on the convergence and are 
the easiest to implement when the grid size is 2 n + 1 where n is any positive in- 
teger. For such a grid it is possible to use n + 1 hierarchical basis levels. Using 
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Figure 3-2: Non-zero Hessian elements for a 9-by-9 image. The plot on the left shows the 
non-zero Hessian elements when the nodal basis is used. The plot on the right shows the 
non-zero elements when a full hierarchical basis is used. 



hierarchical basis functions increases the communication between nodes in the image 
array as shown in Figure 3-2. This increased communication speeds up the diffu- 
sion process considerably. I have noticed a five-fold increase in convergence for the 
photo-topography algorithms. 

The multi-grid methods seek to propagate information over a larger range by 
solving a series of problems of different size. Usually the original problem is formulated 
on grids that decrease in size by a factor of two when going from one layer to the 
next (see Figure 3-3). The solutions on one layer are related to solutions on the layers 
above and below via interpolation or prolongation. The solutions are kept consistent 
with each other via both intra-layer and inter- layer processes (see [Terzopoulos, 1984] 
and [Brandt and Dinar, 1979]). 

Multigrid methods have the potential to be much faster than the hierarchical basis 
functions since most of the computation (and optimization) is done on the smaller 
layers. Multigrid methods are well suited to linear problems (such as surface fitting), 
but may not work for non-linear problems. The problem is that typically 



12 9 (u) ^g(J2u) 
for non-linear functions g(u), and the multi-grid methods require 

12 9 ( u ) =9(12 u ) 



(3.5) 



(3.6) 
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Figure 3-3: Multigrid levels 



to constrain the solutions on the smaller grids so that they don't bias the solution on 
the larger grid. 

One type of multi-grid method that can be used for non-linear problems is the 
coarse-to-fine method. In this method, the problem is solved on coarse layers hrst 
and the solution to each layer provides the initial condition to the next finer layer 
below. This method is significantly faster than just optimizing on the finest grid but 
doesn't produce as much convergence speed-up as the full multi-grid method. 

Like the hierarchical basis methods, the multi-grid methods work best when the 
grid size is 2 n + I . However, since the multi-grid methods define a series of problems 
rather than just choosing a new set of basis functions, the multi-grid methods can be 
implemented easily for all grid sizes. There is some evidence, though, that grid size 
reductions should be near 2 for best convergence rates [Terzopoulos, 1984]. 
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Chapter 4 
Candidate algorithms 



In this chapter I discuss several candidate algorithms for the photo-topography prob- 
lem. Since each of the algorithms estimate the surface depth using a combination of 
shape-from-shading and stereo I call them Depth From Shading and Stereo (DFSS) 
algorithms. Along the way I will discuss the rationale behind each algorithm, its 
strengths and weaknesses, and how it fits in with the fusion techniques discussed so 
far. The performance of these algorithms on a set of synthetic test images is presented 
in Chapter 5. The performance of the z-only algorithm on real images is presented 
in Chapter 9. 

The approach I am taking to fuse the vision algorithms is to develop a single cost 
function that incorporates the problem constraints along with some regularization 
terms to help direct the search path. This approach differs from many researchers' 
attempts at fusion algorithms in that it closely integrates the methods rather than 
building a module-based solution. I believe that this approach will create more robust 
solution methods. 

In effect, I pose a generalized optimization problem which can be solved many 
ways. I have tried various solution methods with this research and have chosen 
direct optimization via the conjugate gradient method as my preferred method. The 
conjugate gradient method, like all direct optimization methods, guarantees reduction 
of the cost function at each step, in contrast to the more widely used Gauss-Seidel 
or Gauss- Jordan methods. The big advantage of the conjugate gradient method over 
other direct optimization methods is that no Hessian needs to be computed or stored. 
Computer vision algorithms have thousands of free variables which would result in a 
Hessian with millions of terms if it were computed. 

I have looked at four basic algorithms for the photo-topography problem. 

• The zpq and z-only algorithms estimate everything in a single global coordi- 
nate system that is defined to be halfway between the two camera positions. I 
refer to these algorithms as centralized. Figure 4-1 shows the flow of a typical 
centralized algorithm as a tree diagram. The current estimate of the surface 
height z is used to project points in the global coordinate system to points in 
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Figure 4-1: Centralized algorithm tree. 
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Figure 4-2: Decentralized algorithm tree. 



each image (via perspective projection). These points won't in general land 
on a pixel center so some type of interpolation (e.g., bilinear or bicubic inter- 
polation) is used to determine the value of the image at the projected points. 
This interpolated image F is then compared to a computed image based on the 
current estimate of the surface. The error is used to update p } q } and ultimately 
z. This type of algorithm closely integrates the constraints from both shading 
and stereo. 

• The dual-z algorithm determines estimates of the surface for each image sepa- 
rately and takes into account the stereo information via a penalty term. I refer 
to this algorithm as decentralized. Figure 4-2 shows the flow of this type of 
algorithm. The estimated height for each image is used to compute the cor- 
responding image estimate. The errors between this estimated image and the 
actual image for each camera are then used to update the associated p } q } and z 
for each image. So far this is the same as for a normal shape-from-shading algo- 
rithm. However, instead of just updating the surface height estimates directly 
from the image errors, the algorithm also takes into account the differences be- 
tween the two surface estimates. This type of algorithm doesn't integrate the 
shading and stereo information as closely as the centralized algorithms. 

• The disparity algorithm is also centralized algorithm but uses estimates of the 
disparity instead of the surface depth as the fundamental variable. Figure 4-3 
shows the flow for this algorithm. The current estimate of the disparity u is used 
to project points in the global coordinate system to points in each image (via 
perspective projection). These points won't in general land on a pixel center 



4.1. ZPQ ALGORITHM. 51 




Figure 4-3: Centralized algorithm tree based on disparity. 

so interpolation is used to determine the brightness value in the image at the 
projected points. This interpolated image F is then compared to a computed 
image based on the current estimate of the surface. The error is used to update 
p, q, and ultimately u. The surface height, z, is computed from u if needed. 
This type of algorithm also closely integrates the constraints from both shading 
and stereo. 

The algorithms differ mostly in the representation of the surface topography. The 
representation chosen for the optimization variables is very important to the per- 
formance of computer vision algorithms. If the right representation is chosen then 
the search space will be smoother and contain less local minima than if the wrong 
representation is chosen. 

4.1 zpq Algorithm. 

This algorithm is based loosely on the Height and Gradient from Shading algorithm 
of Horn [Horn, 1989]. Following Horn's algorithm, surface height z, and the surface 
gradients p } and q, are used as optimization variables. The two photo-topography 
images, camera geometry, and surface reflectance functions are given as inputs to 
the algorithm. The equations that govern this situation are the photo-topography 
constraint equations (2.40) and the integrability constraint equations (2.39). The 
integrability constraint equations are necessary to ensure that p and q are consistent 
with the underlying surface z. The cost function for this algorithm is formed by 
integrating the squared error introduced by the current estimates for p } q } and z, 
together with a penalty functions for departure from integrability and departure from 
smoothness, 

^ J= lJJ)( E{1){x + ?z' y) ~ R{1){P ' q) ) + { E{2) (x-^,y)-^{p,q) > j 
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This cost function allows the normal vector components p and q freedom to try to 
match the input images without requiring integrability. The penalty functions are 
then used to bias p and q toward solutions that are integrable and smooth. Using soft 
constraints like those above provides for more degrees of freedom which, it is hoped, 
will result in fewer local minima. The disadvantages of these extra degrees of freedom 
are a non-exact solution (since z will not be exactly consistent with p and q), and 
slower convergence. 

The cost function above is continuous and must be discretized before it can be 
optimized (solved). The discretization process is an approximation process. The idea 
is to approximate the integral above by some function of a finite (usually small) set of 
variables. Researchers typically use either a series representation of the integrand or 
a sampled version of the integrand. Examples of the hrst representation are Fourier 
series and finite element methods. Examples of the second representation are finite 
difference methods. 

I have used finite difference methods exclusively for the research presented in this 
thesis. The main reason for this is that the images (E^>(x,y) above) are provided in 
digital form. In these images, each pixel represents, in some sense, a weighted average 
of the brightness falling within the sensitive area of the corresponding photosensor. 
Thus each image is actually an array of brightness values. Given this fact it makes 
sense to approximate the values for p } q } and z as arrays of gradient components or 
surface depth, as the case may be. Choosing this digital representation results in a 
problem of the finite difference class. 

Suppose p } q } and z are represented as arrays, as discussed above. Suppose also 
that the continuous derivatives of the underlying functions (such as z) are approxi- 
mated by finite differences, then the cost function can be approximated by the discrete 
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where T> is the discrete domain of the underlying variables in the global coordinate 
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Figure 4-4: Relative sizes of the p, q, and z arrays in comparison to the image arrays for 
the co-grid representation. 

system, M and N are the row and column dimensions of the discrete domain, and e is 
the grid spacing (assumed to the same in both the x and y directions). The F^'(x, y) 
are interpolated from the input images E^'(x,y) using linear interpolation, 



F®(x±£,y) 



E®(x, y) + {x±^-x) [E®(x + l,y)- E^(x, y) 



where 



floor(x ± 



2z> 



(4.3) 



(4.4) 



The floor(x) function returns the greatest integer that is smaller than x. Note that 
Equation 4.3 assumes that x is sampled on a unity-spaced grid. 



4.1.1 Co-grid Implementation 

Two approaches have been taken to implement this cost function. In the hrst ap- 
proach, the p and q arrays are the same size as the image arrays and the z array is 
one pixel larger in both the column and row directions. That is, the p and q functions 
are sampled on the co-grid of the function z (see Figure 4-4). This approach uses 
face-centered surface derivatives that are valid in the center of each 2-by-2 "face" 
formed by the z grid. In this approach, the x- and y- derivatives are approximated by 
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which have low approximation error. The diamond operator o is similar to (but 
different from) two dimensional convolution. 1 In addition, this approach requires 
that z be sampled on the co-grid when computing the terms in the cost function that 
contain z. Based on a bilinear approximation, z can be sampled using 
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The resulting algorithm performs as well as the algorithm to be described next but 
has the disadvantage that the hierarchical basis functions cannot be used to full 
advantage since not all the p } q } and z arrays can be of the size 2 n + 1. Recall that 
the hierarchical basis functions are easiest to implement and provide for optimum 
communication between pixels in the array when the row and column dimensions are 
of size 2 n + I for some n. In addition, multi-grid methods also cannot be used when 
the variable arrays are of different sizes. 

4.1.2 Matched-grid Implementation 

In the second approach, p } q } and z are chosen to all be the same size as the image 
arrays. In such a representation, all the functions are sampled on the same grid (see 
Figure 4-5) which is why this approach is called the matched-grid representation. This 
approach uses vertex- centered surface derivatives that are valid at each vertex of the 



Suppose h(i,j) is a computational molecule (as an array), then the diamond operation is defined 



(hoz)(i,j) = ^2z(i - M + k,j - N + m)h(k,m) 



(4.8) 



where M is the row dimension of h and N is the row dimension of h. The definition differs from 
two-dimensional convolution in that h is not "flipped" . 
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Figure 4-5: Relative sizes of the p, q, and z arrays in comparison to the image arrays for 
the matched-grid representation. 



z grid. Computational molecules that exemplify this approximation are, 
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(4.n; 



which have higher approximation error than their face-centered counterparts. In 
practice, the algorithm based on this approach has nearly the same performance as 
the algorithm based on the face-centered approach. Since p and q are the same size 
as z, some type of approximation must be made at the array edges (boundaries). 
I have chosen to extrapolate the estimates using a bicubic interpolant and use this 
extrapolated version of each estimate in subsequent calculations (see Section A. 2). 
Using the extrapolated estimate is similar to using natural boundary conditions in 
variational calculus methods, and allows for easy computation of the gradient. 

The performance of the matched-grid zpq algorithm on four test image sets is 
presented in Chapter 5. This algorithm gets stuch in a local minimum on the hard 
crater images (as do all of the algorithms to be presented), but does reasonable well 
on the other test images (except for the mountain images). As a whole, this algorithm 
produces surfaces that are too smooth and tahes longer to converge than the other 
algorithms. 
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Like all the algorithms to be presented, this algorithm is run using a exponentially 
changing smoothness parameter A, and integrability parameter // (if it exists). A 
large initial smoothness parameter helps to make the search space more bowl-like 
at the expense of meeting the photo-topography constraints. When the smoothness 
parameter is small, the photo-topography constraints are met at the expense of a 
possibly bumpy surface. A small initial integrability parameter allows the surface 
gradient components p and q to wander far from integrability in order to create image 
estimates that match the input images. A large integrability parameter, requires that 
these image estimates be created so that p and q are nearly integrable. 

While it is not required that the optimization parameters be ramped in this fash- 
ion, I have found that doing so allows for faster convergence over optimizing with the 
parameters fixed at their final values. 

4.2 2>only Algorithm 

This next cost function is similar to the zpq cost function except that it uses hard 
integrability constraints instead of soft constraints. With hard constraints we are 
guaranteed that any solutions obtained will be feasible. The trade-off is that the 
algorithm will have less degrees of freedom to work with and may be more susceptible 
to local minima. 

For this cost function, the depth map z is the only optimization variable thus the 
resulting algorithm is called the z-only algorithm. The surface gradient components 
p and q are computed directly from the depth map. As before, the photo-topography 
images, camera geometry, and surface reflectance functions are inputs to the algo- 
rithm. The constraint equations again come from the photo-topography constraints 
and integrability constraints. The cost function is formed by integrating the squared 
photo-topography error introduced by the current estimate for z, together with a 
penalty function for departure from smoothness, 

The penalty function is mainly used to guide the solution towards the minimum. 
In practice, the smoothness weighting parameter, A, is slowly reduced toward zero as 
the algorithm converges. 

+ \[z 2 xx + 2z 2 xy + z 2 yy ]}dxdy. (4.12) 

The smoothness term is based on what is called the second variation. It is equivalent 
to p x + p 2 + q 2 + q 2 when z x sa fp/z and z y sa fq/z where z is the nominal depth. 2 



2 The approximations are valid when the field of view is small, the image is centered around the 
camera's principal point, and the depth of field relative to the nominal depth is small. 
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Figure 4-6: Relative size of the z array in comparison to the image arrays for the expanded 
boundary implementation. 

Using an array for z, the discrete approximation of this cost function is, 
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+ \FW(x-£,y)-I*Q(p,q ) 

+ \[z 2 xx + 2z 2 xy + z 2 yy }}. 



(4.13) 



The F^'(x± |^, y) are interpolated from the input images E^'(x, y) as in the previous 
algorithm. 



4.2.1 Expanded Boundary Implementation 

Two approaches have been taken to implement this cost function. In the hrst ap- 
proach, the z array is 1 pixel larger than the image array in all directions (see Fig- 
ure 4-6). The normal vector components are estimated on the inner grid which is the 
same size as the image arrays. In fact, the outer depth estimates are only used when 
computing p and q. This approach is very easy to implement but cannot be used with 
a multi-grid scheme since the z array is a different size than the image arrays. The 
approach can, however, take full advantage of hierarchical basis functions since z can 
be of the size 2 n + 1. The implementation uses the vertex-centered derivatives dis- 
cussed in the previous algorithm. The second-order derivatives are approximated by 
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the following computational molecules which are also vertex- centered approximations, 
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The algorithm based on this approach has nearly the same performance as the algo- 
rithm based on the approach to be presented next. 



4.2.2 Matched Grid Implementation 

In the second approach, z is chosen to be the same size as the image arrays thus 
facilitating a multi-grid implementation. The surface normal components, p and q, 
are computed on the boundary using the boundary extrapolation technique discussed 
for the zpq algorithm. This approach also uses vertex-centered approximations for 
the derivatives. 

The performance of the matched-grid z-only algorithm on the four test images 
is presented in Chapter 5. This algorithm also gets stuck in the local minimum for 
the hard crater images. The performance on the easy crater images and the others 
is very good however. In fact, this algorithm has the best performance of the all the 
algorithms to be presented. An important thing to note in the performance figures 
is that the algorithm converges to a reasonable surface quite quickly, usually in the 
hrst 100 iterations or so. The performance figures are shown for 1200 iterations so 
that it is possible to see the convergence characteristics in the long term. 
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4.3 Dual-^ Algorithm 

This cost function uses a totally different variable representation. In this case the 
depth map is represented by two arrays, z^- 1 ' and z^ 2 ' which are estimated based 
on the images E^ 1 ' and E^ 2 > respectively. This algorithm is referred to as the dual-z 
algorithm. The two depth arrays are kept consistent with one another using a penalty 
function that enforces the stereo constraint equations, 



Z W {x MyW) =Z W {x W-J£,yi 



z( 



(4.18) 



where (x^ 1 ' , y^ 1 ') are the coordinates of points in the hrst image and (x^ 2 \ y^ 2 ') are the 
coordinates of points in the second image. These equations require that the depth 
estimate from one image be the same as the depth estimate at the corresponding 
projected point 3 in the other image. Each depth map is a constrained shape-from- 
shading solution for the corresponding image. This algorithm has a totally different 
search space than any of the previous algorithms. 

For this cost function, the depth maps z^> are the optimization variables and 
the surface gradients p^' and q^' are computed directly from the z^> . The photo- 
topography images, camera geometry, and the surface reflectance functions are inputs. 
The constraints come from the photo-topography constraints, the integrability con- 
straints, and the stereo constraints. The cost function is formed by integrating the 
squared shading error (in the spirit of shape-from- shading) introduced by the current 
estimates of z^\ together with penalty functions for departure from smoothness, and 
stereo error. 



mm 



J=\jj{ (£ ( V x) , y (1) ) - R {1) ( P {1 \q {1) )) 2 + X(V 2 z^) 2 

+ // (.W(xW, j/W) - z^(x^ - ^, yW)\ 2 1 dxW dyW 

+ l -jj{ (£ (2 V 2) , y (2) ) - i? (2 V 2 W 2) )) 2 + A(vV 2 )) 2 

+ // (z^(x^ + ||, yW) - z^(x^\ yW)) ' 1 dx™ dy^ (4.19) 



where 



V 2 z = z xx + z yy . (4.20) 

The smoothness penalty function for this cost function is based on the squared Lapla- 



5 Using perspective projection. 
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cian of z^>. The squared Laplacian has roughly the same characteristics as the second 
variation [Grimson, 1979], but is much quicker to compute (since it only involves one 
convolution). 

The cost function contains two shape-from-shading cost functions and the stereo 
penalty function. This cost function is not as closely-coupled as the other two cost 
functions I have described so far. In particular the stereo constraints are added on 
like an addendum to the shape-from-shading cost function. The stereo constraints, in 
effect, pass their constraints up to the shape from shading parts of the algorithm. I 
think of this algorithm as trying to fit together two rubber mountains. The mountains 
can be moved back and forth, and molded to find the best fit. Since the algorithm 
is not as integrated, we would expect that images that don't contain strong shading 
information will be difficult for this algorithm to solve. 

Using z^- 1 ' and z^ 2 > arrays as before, the cost function can be discretized to become, 
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The stereo penalty functions include terms involving z^ 1 ' and z^ 2 ' which are inter- 
polated versions of z^' and z^ 2 '. Bicubic interpolation is used [Keys, 1981]. Using 
this notation, z^'(x,y) is the value of z^> evaluated at the point (x,y) using bicubic 
interpolation. Vertex-centered derivative approximations were used with this cost 
function, and the Laplacian is approximated using the computational molecule, 
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Both the expanded boundary approach and the matched-grid approach discussed 
in the previous section were taken to implement this cost function. Several variations 
of the basic implementation were tried with this cost function. For example, in 
addition to the smoothness term based on the squared Laplacian, algorithms were 
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implemented where the smoothness term was based on the second variation, 

+ ^W S {■ ■ ■ + A [(4 2 2) 2 + 2(4?) 2 + (4?) 2 ] • • • +} -(4-23) 

Another alternate implementation used face-centered derivatives instead of vertex- 
centered derivatives. All of these algorithms have performance very similar to each 
other, and they suggest that as long as a reasonable derivative approximation or 
smoothness criteria is used, it doesn't really matter which one is chosen. Any reason- 
able approximation or smoothness criteria is as good as any other. In light of this, it 
makes sense to choose the approximations or smoothness criteria that are the easiest 
and the fasted to compute, or that have the best numerical properties. 

The performance of this algorithm on four test images is presented in Chapter 5. 
Like the previous algorithm, this algorithm also has problems with the hard crater 
images but does well on the other images. Unfortunately, the results of this algorithm 
are two depth maps which represent the depths as seen from the right and left cameras. 
These maps will not, in general, register exactly with one another. To create a 
single depth map from these maps requires some type of averaging. This is a major 
disadvantage. 

To address this disadvantage, an algorithm was implemented that had three depth 
maps z^ 1 ' , z( 2 \ and z, where z is the depth map as seen from the global coordinate 
system. In all other ways the algorithm was the same as the dual-z algorithm. The 
performance of the algorithm was very similar to the performance of the above al- 
gorithm yet it created a central depth map. It is not shown here since it was even 
slower than the dual-z algorithm. 

The combination of two depth maps and need for bicubic interpolation, caused 
the dual-z algorithm to be the slowest to compute. The figures in Chapter 5 show 
the number of function evaluations and not the computation time. In practice this 
algorithm was 3-5 times slower than the z-only algorithm. 

4.4 Disparity Map Algorithm 

Another popular representation for depth that shows up often in the vision literature 
is the disparity map. The disparity is the relative offset of a point in one image with 
respect to the corresponding point in the other image. Using the notation we have 
used so far, the disparity can be defined as, 

u = g, (4,24) 
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where u is the disparity. Using this definition of disparity, the surface gradient com- 
ponents are computed using, 

P= — , 



q= J y . 

yu y - u 

The photo-topography constraint equations can also be re-written using disparity to 
become, 

E^(x + u,y) = R^(p,q), 

EW(x-u,y) = R?)(p, q ). (426) 

Suppose we adopt the centralized approach for this cost function and have the dis- 
parity map u be the optimization variables. The surface gradients p and q are com- 
puted directly from u. The photo-topography images, camera geometry, and surface 
reflectance functions are inputs. The cost function for this case is formed by integrat- 
ing the squared photo-topography error (based on disparity) together with a penalty 
function for departure from smoothness, 

nAnJ = ^JJ{(E( 1 \x + U ,y)-RW(p,q))\(E( 2 \x-u,y)-RW(p,q)] 

+ \[ul x + 2ul y + u 2 yy ]}dxdy. (4.27) 

where the surface gradient components are given by Equation 4.25. The smoothness 
term is based on the second variation of u. The main difference between this cost 
function and the others I have presented is the smoothness term (it is now based on 
disparity instead of depth). The numerical properties will also be somewhat different 
since the hierarchical basis functions will be over disparity instead of depth. 

Disparity is typically used in vision algorithms for two reasons. First the disparity- 
based algorithms usually have fewer division operations which can improve perfor- 
mance, and second the disparity has a much smaller dynamic range which can con- 
tribute to better numerical properties. Unfortunately, since the photo-topography 
problem requires the use of perspective projection, the hrst reason doesn't hold for 
my algorithms. The culprits are the surface gradient equations (Equations 4.25). As 
for the second reason, the disparity-based DFSS algorithm has pretty much the same 
numerical properties as the other DFSS algorithms. However, using disparity, instead 
of depth, changes the solution space the optimization is performed in. 

Representing u as an array, the cost function can be discretized to obtain, 

m « in J = 2mW ^ {{ F(1) ( x + u ,y)~ R {1) (p,q)) 2 + {f {2) (x -u, y )- R {2 \ P , q ) " 

x,yEV 



+ A [t4 + 2<, + «;J } (4.28) 
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The F^'(x ± u,y) are linearly interpolated from the input images E^'(x,y) using 
disparity, 

F (i) {x±u,y) = E (l \x,y) + {x±u-x) [E^(x + l,y) - E (l \x, y)\ (4.29) 

where 

x = floor(x±u). (4.30) 

Only the matched-grid approach was used to implement this cost function. In 
this case, u is chosen to be the same size as the image arrays. The derivatives 
are approximated using vertex-centered computational molecules as for the previous 
algorithms. 

The performance of this algorithm on four test images is presented in Chapter 5. 
This algorithm, like all the previous algorithms, gets stuck in a local minimum for 
the hard crater images. The algorithm, however, performs well on the other images. 
An interesting observation can be made about the performance of this algorithm on 
the mountain images. The algorithm over-estimates the mountain height while the 
z-only algorithm under-estimates it. This is clear indication that the search space for 
the disparity-based algorithm is fundamentally different from the search space for the 
z-only algorithm. Unfortunately, the solution space is not different enough to avoid 
falling into the local minimum of the hard crater image. 
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Chapter 5 
Test Results 



It is important when developing a new algorithm to be able to test its performance. 
For computer vision algorithms, that means that it must be possible to compare the 
estimated surface with the actual surface. The only way to do this with complete 
confidence in the results is to create synthetic images from a known surface topogra- 
phy, run the algorithm on these images, and compare the estimate with the known 
topography. Only after it is shown that the algorithm works on the synthetic im- 
ages, can we be confident that the surface estimates for actual planet images will be 
accurate. 

I have created four image pairs as test images based on three surface topologies. 

• Easy Crater Images. This set of images is based on a crater on a flat plane. 
The light sources are oblique which results in images that prove to be easy for 
the DFSS algorithms to interpret. 

• Hill Images. This set of images is based on a fractally-generated set of rolling 
hills. The light sources are oblique. These images are interpreted correctly for 
most of the DFSS algorithms but require more iterations of each algorithm since 
they are more complicated than the crater images. 

• Mountain Images. This set of images is based on a fractally-generated moun- 
tainous terrain. The light sources are oblique and the camera baseline is much 
smaller. This set of images poses a challenge to the DFSS algorithms due to 
the steep terrain and reduced baseline. 

• Hard Crater Images. This set of images is also based on a crater on a flat 
plane. The light sources are almost directly behind the camera resulting in a 
set images that prove to be difficult for the DFSS algorithms. 

The calibration parameters for the test images are summarized in Table 5-1. The 
table lists values for the baseline distance 6, camera focal distance /, nominal depth 
z , and light source vector components p s , q s , p s , q s . The focal distance and 
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Table 5-1: Camera geometry. 

nominal depth are negative to be consistent with a right handed coordinate system. 
The baseline distance b and depth z have the same units (say miles), while the camera 
local distance, pixel spacing and the light source components are based on camera 
units (say millimeters). The table also lists values lor relative surface height Az/z 
and disparity Am. These two parameters indicate the difficulty ol the problem. Small 
values ol either indicate a hard problem. Note the units lor Am are in pixels lor a 
65-by-65 problem, and / is computed so that the pixel spacing is 1 (millimeter). 

All of the test images are noise-free so that the best performance of each algorithm 
can be tested. The performance of the z-only algorithm on images with noise in 
presented in Chapter 7. The performance figures for each test case include mesh plots 
of the estimated surface at several points during optimization as well as estimated 
images based on the final surface. The mesh plots are of a 33-by-33 smoothed (to 
avoid aliasing) and subsampled version of the 65-by-65 estimated surface. The lower 
resolution mesh plots are used to avoid printing problems. 

The performance figures also show the number of function evaluations (iterations) 
computed by each algorithm during convergence. The number of conjugate gradient 
updates taken is between one-half to one-third of the number of function evaluations 
since 2-3 function evaluations are needed to perform each update step. All of the 
algorithms use hierarchical basis functions to enhance convergence. 

5.1 Easy Crater Images 

The crater on a flat plane is very simple surface and thus serves as a good test 
surface. 1 The camera and light source geometry as well as the true surface, light 
source contour plots and resulting images are shown in Figure 5-1. As shown in the 
figure, the baseline distance between the cameras is about half the distance to the 
surface and the light source (i.e., the sun) positions for the two images differ greatly 
from each other. As a result, this set of images has strong shading information which 
is very easy for the DFSS algorithms to take advantage of. The effects of the lighting 
can be seen in the contour plots of the reflectance map. The reflectance contours 
are separated enough in gradient space (p-q space) so that brightness values from 



information on how to generate this surface can be found in Appendix C. 
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Figure 5-1: Test images of crater on flat plane (easy case). Shown are the camera geometry 
as projected into the xz- and y^-planes, the true surface as a mesh plot, reflectance function 
contours for the two light source positions, and the left and right synthetic noise-free images. 



68 CHAPTERS. TEST RESULTS 

one image can easily constraint the possible set of feasible gradient directions. 2 The 
images are 65-by-65 pixels in size and are noise free. 

The Figures 5-2-5-5 show the results of applying the various DFSS algorithms 
to this test case. All four algorithms can correctly interpret these images but the 
z-only algorithm performs best. The figures show that the estimated images in this 
case very closely resemble the true images and the surface estimates are very good 
except in the lower left corner of the surface (see the mesh plots). With the lighting 
conditions chosen for these test images, this small anomaly has little effect on the 
estimated images. It is interesting to note that the anomaly shows up for all four 
algorithms. 

The most interesting thing about this test case is the rate of convergence that is 
obtained with the z-only, dual-z, and disparity algorithms. In those cases, a pretty 
good estimate is obtained after only 50 function evaluations! This represents only 
about 20 updates since more than one function evaluation is necessary per update 
when using conjugate gradient optimization. 



2 That is, there are regions in the gradient space where a given brightness value from one image 
strongly constrains the gradient direction there. 
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Figure 5-2: Performance of the zpq algorithm on the easy crater images. The up- 
per-left graph shows the cost function (solid line) and RMS error of the estimated surface 
(dashed-line) plotted against the number of function evaluations. The upper-right graph 
shows the history of the A (solid line) and fj, (dashed-line) cost function parameters. Also 
shown are images based on the estimated surface and mesh plots of the surface at different 
points during optimization. 
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Figure 5-3: Performance of the z-only algorithm on the easy crater images. The up- 
per-left graph shows the cost function (solid line) and RMS error of the estimated surface 
(dashed-line) plotted against the number of function evaluations. The upper-right graph 
shows the history of the A cost function parameter. Also shown are images based on the 
estimated surface and mesh plots of the surface at different points during optimization. 
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Figure 5-4: Performance of the dual- z algorithm on the easy crater images. The up- 
per-left graph shows the cost function (solid line) and RMS error of the estimated surface 
(dashed-line) plotted against the number of function evaluations. The upper-right graph 
shows the history of the A (solid line) and fj, (dashed-line) cost function parameters. Also 
shown are images based on the estimated surfaces and mesh plots of the surfaces at different 
points during optimization. 
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Figure 5-5: Performance of the disparity-based algorithm on the easy crater images. The 
upper-left graph shows the cost function (solid line) and RMS error of the estimated surface 
(dashed-line) plotted against the number of function evaluations. The upper-right graph 
shows the history of the A cost function parameter. Also shown are images based on the 
estimated surface and mesh plots of the surface at different points during optimization. 
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5.2 Hill Images 

The second set of images is of an undulating surface similar to eroded hills and was 
generated using a fractal technique. 3 Figure 5-6 shows the camera geometry, true 
surface, reflectance contours, and test images for this surface. This set of images is 
more representative of the type of terrain the DFSS algorithms are likely to encounter. 
As shown in the figure, the left camera is directly over the surface and the right camera 
views the surface obliquely. The light sources are separated, as for the easy crater 
image, resulting in images with strong shading information. I find the relatively low 
65-by-65 resolution of these images hard to interpret visually (unlike the crater images 
which are easy to interpret). The DFSS algorithms, however, perform well with these 
images. 

The Figures 5-7-5-10 show the results of applying the various DFSS algorithms 
to this test case. All four algorithms work reasonably well with these images, but 
the z-only algorithm performs best. Due to the complexity of the surface, this set of 
images requires more iterations of each algorithm to obtain a satisfactory estimate of 
the surface than the easy crater images. The figures show that the z-only algorithm 
performs best, it even obtains a very good estimate after 50 iterations! The con- 
vergence is a little slower for the dual-z and disparity algorithms; they obtain good 
estimates after about 200 iterations. The zpq algorithm is the slowest and requires 
about 600 iterations to obtain a satisfactory solution. The estimated images in all 
cases match very well the true images shown in Figure 5-6. 



information on how to generate this surface can be found in Appendix C. 
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Figure 5-6: Test images of hill. Shown are the camera geometry as projected into the xz- 
and y^-planes, the true surface as a mesh plot, reflectance function contours for the two 
light source positions, and the left and right synthetic noise-free images. 



5.2. HILL IMAGES 



75 



10 



10 



s 10- 



10 



Convergence History 



500 1000 

Iteration 



Lambda and Mu Parameter History 



10 



1500 




500 1000 

Iteration 



1500 




Estimated images at iteration 1401 



Estimated surface at iteration 50 



Estimated surface at iteration 201 





Estimated surface at iteration 600 



Estimated surface at iteration 1401 





Figure 5-7: Performance of the zpq algorithm on the hill images. The upper-left graph shows 
the cost function (solid line) and RMS error of the estimated surface (dashed-line) plotted 
against the number of function evaluations. The upper-right graph shows the history of the 
A (solid line) and fi (dashed-line) cost function parameters. Also shown are images based on 
the estimated surface and mesh plots of the surface at different points during optimization. 
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Figure 5-8: Performance of the z-only algorithm on the hill images. The upper-left graph 
shows the cost function (solid line) and RMS error of the estimated surface (dashed-line) 
plotted against the number of function evaluations. The upper-right graph shows the history 
of the A cost function parameter. Also shown are images based on the estimated surface 
and mesh plots of the surface at different points during optimization. 



5.2. HILL IMAGES 



77 



Convergence History 




Lambda & Mu Parameter History 



500 1000 1500 

Iteration 




500 1000 1500 

Iteration 




Estimated images at iteration 1201 



Estimated surfaces at iteration 5 1 



Estimated surfaces at iteration 201 




Estimated surfaces at iteration 600 



Estimated surfaces at iteration 1 20 1 



Hi 




Figure 5-9: Performance of the dual- z algorithm on the hill images. The upper-left graph 
shows the cost function (solid line) and RMS error of the estimated surface (dashed-line) 
plotted against the number of function evaluations. The upper-right graph shows the history 
of the A (solid line) and fj, (dashed-line) cost function parameters. Also shown are images 
based on the estimated surfaces and mesh plots of the surfaces at different points during 
optimization. 
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Figure 5-10: Performance of the disparity-based algorithm on the hill images. The up- 
per-left graph shows the cost function (solid line) and RMS error of the estimated surface 
(dashed-line) plotted against the number of function evaluations. The upper-right graph 
shows the history of the A cost function parameter. Also shown are images based on the 
estimated surface and mesh plots of the surface at different points during optimization. 
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5.3 Mountain Images 

The third set of images is of a highly mountainous surface and was created to show 
the performance of the algorithms on steep terrain. The steep terrain requires that 
the baseline distance from the cameras be shortened so that all the surface points 
are visible from both cameras (see Figure 5-11). Thus this set of images can be 
used to test the performance of the algorithms when the stereo baseline is small. 
As shown in the figure, the light source positions for this set of images are widely 
separated and generate deep shadows on this steep terrain. The reflectance maps 
are flat within a shadow so no helpful gradient is available to the algorithms. In 
addition, knowledge that a particular pixel is in shadow only constrains the set of 
possible gradient directions to a sub-plane of gradient space. Thus, within a shadow 
region, much more influence is given to the brightness values from the other image. 
This combination of effects results in slower convergence. 

The Figures 5-12-5-15 show the results of applying the various DFSS algorithms to 
this test case. The z-only, dual-z, and disparity algorithms work reasonably well with 
these images, but the z-only algorithm performs best. The zqq algorithm has trouble 
with this test case. Due to the complexity of this test case, the algorithms require 
many more iterations to achieve a satisfactory solution. The z-only and disparity 
algorithms are able to obtain good solutions after about 600 iterations, while 1400 
iterations are not enough for the zpq. The dual-z algorithm obtains a reasonable 
estimate at about 1200 iterations. 
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Figure 5-11: Test images of mountain. Shown are the camera geometry as projected into 
the xz- and y^-planes, the true surface as a mesh plot, reflectance function contours for the 
two light source positions, and the left and right synthetic noise-free images. 
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Figure 5-12: Performance of the zpq algorithm on the mountain images. The up- 
per-left graph shows the cost function (solid line) and RMS error of the estimated surface 
(dashed-line) plotted against the number of function evaluations. The upper-right graph 
shows the history of the A (solid line) and fj, (dashed-line) cost function parameters. Also 
shown are images based on the estimated surface and mesh plots of the surface at different 
points during optimization. 
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Figure 5-13: Performance of the z-only algorithm on the mountain images. The up- 
per-left graph shows the cost function (solid line) and RMS error of the estimated surface 
(dashed-line) plotted against the number of function evaluations. The upper-right graph 
shows the history of the A cost function parameter. Also shown are images based on the 
estimated surface and mesh plots of the surface at different points during optimization. 
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Figure 5-14: Performance of the dual-,? algorithm on the mountain images. The up- 
per-left graph shows the cost function (solid line) and RMS error of the estimated surface 
(dashed-line) plotted against the number of function evaluations. The upper-right graph 
shows the history of the A (solid line) and fj, (dashed-line) cost function parameters. Also 
shown are images based on the estimated surfaces and mesh plots of the surfaces at different 
points during optimization. 
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Figure 5-15: Performance of the disparity-based algorithm on the mountain images. The 
upper-left graph shows the cost function (solid line) and RMS error of the estimated surface 
(dashed-line) plotted against the number of function evaluations. The upper-right graph 
shows the history of the A cost function parameter. Also shown are images based on the 
estimated surface and mesh plots of the surface at different points during optimization. 
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5.4 Hard Crater Images 

Like the easy crater images presented earlier, this set of images is also based on the 
crater on a flat plane. However in this case, the light source positions are nearly the 
same and are almost directly behind the cameras. The resulting images are very bright 
and look very much the same except for slight differences in shading (see Figure 5- 
16). The reflectance maps for this case are nearly co-incident which contributes to 
a surface orientation ambiguity. This set of images represents a worst case for the 
DFSS algorithms since the shading information in the images is weak and the range 
of brightness values is small. However, they also have strong stereo correspondence 
information, which unfortunately is not readily utilized by the DFSS algorithms. The 
test images are 65-by-65 pixels in size and are noise free. 

The Figures 5-17-5-20 show the results of applying the various DFSS algorithms 
to this test case. All four algorithms get stuck in a local minimum. The figures 
show that the estimated images resemble the true images of Figure 5-16 even though 
the estimated surface doesn't match the actual surface. The problem is that the 
algorithms incorrectly interpret the surface as concave when it is actually convex (see 
also the shape-from-shading results in Chapter 6). Even though the stereo information 
in these test images can be used to correctly determine the orientation of the surface, 
the DFSS algorithms rely very heavily on the shading information. 

The surface orientation ambiguity is a result of having both light sources directly 
behind the cameras. All the test images that I have tried that have this light source 
geometry cause the algorithms to fail. 

The results shown in the figures are based on algorithms that employ the hier- 
archical basis functions. When I run the algorithms on the hard crater test images 
without using the hierarchical basis functions, the algorithms can correctly interpret 
the images. The convergence is significantly slower, however. Based on several runs I 
have made, it appears that the explanation for this behavior is that the hierarchical 
basis functions lead the algorithm down a particular path that ends up in the local 
minima. When the nodal basis is used, the algorithm goes down a different path that 
is able to bypass the local minimum (just barely). 

This problem seems to be sensitive to initial condition. In fact, when I start the 
algorithms from an initial condition close to the true surface, they do converge to the 
true solution as expected. 4 

This test case has been the impetus for developing a plethora of algorithms, only 
a few of which are presented in this thesis. I had hoped that developing an algorithm 
with a totally different search space would result in one that could solve this problem. 
Unfortunately, all the algorithms I have developed have gotten stuck in the same local 
minima shown in the figures. The other test cases show that the algorithms perform 
very well when different lighting is used where at least one lighting condition is oblique. 



4 Of course, they only achieve the exact solution when the smoothness parameters are set to zero. 
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Figure 5-16: Test images of crater on flat plane (hard case). Shown are the camera geometry 
as projected into the xz- and y^-planes, the true surface as a mesh plot, reflectance function 
contours for the two light source positions, and the left and right synthetic noise-free images. 
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Figure 5-17: Performance of the zpq algorithm on the hard crater images. The up- 
per-left graph shows the cost function (solid line) and RMS error of the estimated surface 
(dashed-line) plotted against the number of function evaluations. The upper-right graph 
shows the history of the A (solid line) and fj, (dashed-line) cost function parameters. Also 
shown are images based on the estimated surface and mesh plots of the surface at different 
points during optimization. 
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Figure 5-18: Performance of the ^-only algorithm on the hard crater images. The up- 
per-left graph shows the cost function (solid line) and RMS error of the estimated surface 
(dashed-line) plotted against the number of function evaluations. The upper-right graph 
shows the history of the A cost function parameter. Also shown are images based on the 
estimated surface and mesh plots of the surface at different points during optimization. 
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Figure 5-19: Performance of the dual-2; algorithm on the hard crater images. The up- 
per-left graph shows the cost function (solid line) and RMS error of the estimated surface 
(dashed-line) plotted against the number of function evaluations. The upper-right graph 
shows the history of the A (solid line) and fj, (dashed-line) cost function parameters. Also 
shown are images based on the estimated surfaces and mesh plots of the surfaces at different 
points during optimization. 
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Figure 5-20: Performance of the disparity-based algorithm on the hard crater images. The 
upper-left graph shows the cost function (solid line) and RMS error of the estimated surface 
(dashed-line) plotted against the number of function evaluations. The upper-right graph 
shows the history of the A cost function parameter. Also shown are images based on the 
estimated surface and mesh plots of the surface at different points during optimization. 
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Absolute Error 
Algorithm Easy Crater Hill Mountain Hard Crater 



dual- 2; 


1.2991 


0.2274 


0.8774 


2.4241 


zpq 


0.2130 


0.0247 


2.7494 


1.5933 


disparity 


0.1668 


0.0251 


0.8383 


0.9241 


z-only 


0.1728 


0.0247 


0.7799 


1.0716 


Algorithm 


Relative Error 
Easy Crater Hill Mountain 


Hard Crater 


dual- 2; 


0.3614 


0.0968 


0.8436 


1.0078 


zpq 


0.1678 


0.0205 


2.5715 


1.5379 


disparity 


0.1663 


0.0176 


0.7193 


0.9240 


z-only 


0.1721 


0.0175 


0.7177 


1.0078 



Table 5-2: Absolute and relative surface estimation error by algorithm and test case (1200 
iterations). 



5.5 Summary 

A summary of running the algorithms on the test images is shown in Table 5-2 and 
Figure 5-21. The table shows the relative and absolute error between the true and 
estimated surface at the last iteration. The absolute error is computed using the 
formula 



J. 
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while the relative error is computed using 



J, 
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(5.2) 
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where z* is the true surface height and z is the average of z. The figure shows 
the absolute error normalized by the surface depth change of the true surface (i.e., 
{z max — z m i n )). The figure can be used to interpret the error as a fraction of the total 
surface depth change. All the algorithms were run with 6 hierarchical basis levels 
which is one less than the largest number of levels that can be used with 65-by-65 
images. Using this number of hierarchical bases sped up the convergence of each 
algorithm by a factor of 3-5. 

It is clear that all of the algorithms had problems with the hard crater images. 
Each one got caught in the local minimum. The reason for this can be seen in the 
reflectance contours of Figure 5-16. For a given brightness level (i.e., along one of 
the contours), there are two viable solutions with different surface orientations. The 
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Figure 5-21: Absolute surface estimation error scaled by the true surface depth change. The 
depth change is shown in parentheses. 



local minimum has a orientation in the "dipped" region that is a viable but incorrect 
orientation, and the stereo information is not strong enough to pull it out of the local 
minimum. 

The hard crater images have been the driving force behind much of this research 
and are one of the reasons that so many alternate algorithms were implemented. 
As each algorithm was created and implemented it was hoped that the new search 
space would allow the correct interpretation of the hard crater images. Alas, every 
algorithm I have tried has been caught in the local minimum. 

The only way the algorithms can correctly interpret the hard crater images is to 
perform the optimization with no hierarchical basis levels. The convergence is very 
slow but the solution is correct. I hesitate to recommend that all runs be performed 
using the nodal basis, however, since the performance hit is huge and there is no 
guarantee that the true solution will be found. Consider the fact that the test images 
are small and that more normal size images are of the order 513-by-513. One function 
evaluation with these images will take 64 times as long to compute as the test images. 
In addition, since all of the algorithms are diffusion-type algorithms, we could expect 
an 8 times increase in the number of function evaluations to obtain convergence if the 
nodal basis is used. As you can see, the ramifications of not using the hierarchical 
bases are huge. 

Fortunately, the fact that the algorithms are caught in a local minimum is easy 
to see for the hard crater images, since the estimated images that are created as 
a by-product of each function evaluation show streaks where the photo-topography 
constraints are not met (these streaks are difficult to see in the halftoned reproductions 
of the gray level images). Considering this fact, I recommend that the largest number 
of hierarchical basis levels be used to speed convergence and the results checked by 
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looking at the estimated images. 

Clearly, the algorithm that performed best is the z-only algorithm. Not only is 
it the quickest algorithm to compute (i.e., one iteration of the z-only algorithm com- 
putes faster that one iteration of any of the other algorithms), 5 but it also has the 
best convergence rate. This algorithm and the next best algorithm, the disparity 
algorithm, were created by directly implementing the problem constraints as a cost 
function. There might be a lesson here that the best algorithms are formed by in- 
corporating the constraint equations into the cost function in the most direct and 
straightforward way. 



5 The evaluation time for the algorithms is approximately proportional to the number of free 
variables; N 2 for z-only and disparity, ?>N 2 for zpq, and 2N 2 for dual-z. 
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Chapter 6 

Comparison with 
Shape- from- Shading 



In this chapter I show how a simple shape- from-shading algorithm performs on my test 
images. By restricting the z-only algorithm to work on a single image and assuming 
orthographic projection I can create an algorithm similar to Szeliski [Szeliski, 1991] 
and Leclerc and Bobick [Leclerc and Bobick, 1991]. The resulting algorithm is based 
on the cost function, 

min J=\JJ {(E(x, y) - R(p, q)f + \{V 2 z) 2 } dx dy (6.1) 

where 

P = Z X1 , 

F ' (6.2) 

q = z y . 

The shape- from-shading algorithm is based on z but cannot estimate the true depth; 
the algorithm can only estimate the shape. The depth bias is unobservable to this 
algorithm. 

Following the development of the z-only algorithm, the discretized version of this 
cost function is 

mm J = ^^ £ {(E(x, y) - R(p, q)f + X(V 2 zf} . (6.3) 

x,y(zV 

Note that this cost function directly uses the values in E(x } y) without having to 
interpolate as was required when forming F(x } y) for the z-only algorithm. 

I show the performance of the shape-from-shading algorithm on the left image of 
four test image sets in Figures 6-4-6-3. The figures show the cost function history 
and optimization parameters as a function of the number of function evaluations. 
Also shown is the estimated surface shape at various stages during the convergence. 
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Figure 6-1: Performance of the shape-from-shading algorithm on the left easy crater image. 
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Figure 6-2: Performance of the shape-from-shading algorithm on the left hill image. 
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Figure 6-3: Performance of the shape-from-shading algorithm on the left mountain image. 



99 



10 



Convergence History 



10" 



10 



Lambda Parameter History 



500 1000 

Iteration 



1500 




500 1000 1500 

Iteration 




Estimated image at iteration 1 200 



Estimated surface at iteration 50 



Estimated surface at iteration 200 





Estimated surlace at iteration 600 



Estimated surface at iteration 1200 





Figure 6-4: Performance of the shape-from-shading algorithm on the left hard crater image. 
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Image set Rel. SFS Error Rel. z-only Error 



Easy Crater 


0.6242 


0.1721 


Hill 


0.2083 


0.0175 


Mountain 


1.6343 


0.7177 


Hard Crater 


1.3532 


1.0715 



Table 6-1: Estimation error of SFS algorithm on the test images compared to z-only results 
(1200 iterations). 



The figures show that the shape-from- shading algorithm can determine the shape 
of the surface, but with some error. I have summarized the final error between the 
estimated shape and true shape in Table 6-1. The estimated surface and true surface 
offset so that their average heights coincide before computing the errors shown in the 
table. The error was computed using Equation 5.2. 

The estimation error when using shape-from-shading alone is greater than the 
estimation error when using the z-only algorithm. While this is not surprising since 
the z-only algorithm has two images to work with instead of one, it is nice to see that 
significant performance gains are possible when using a fused algorithm. 



Chapter 7 
Error Analysis 



It is equally important when developing algorithms to investigate their robustness in 
the face of errors. In this chapter I present the results of running the z-only algorithm 
on hill images that contain errors. I have introduced errors of three types: 

1. measurement errors, 

2. geometry errors, 

3. reflectance errors. 

Each error is introduced separately to simplify the analysis. In each case, a figure is 
presented that shows the relative and absolute error between the true and estimated 
surface, along with the estimated surface, convergence history and estimated images 
for the worse case. All the results below are the error after 600 function evaluations 
(about 300 updates). The initial condition for each case was a flat plane at the 
nominal depth. For reference, the relative and absolute errors after 600 function 
evaluations on noise-free images are 0.0175 and 0.0247 respectively. 

7.1 Measurement Errors 

To determine the performance of the z-only algorithm on noisy images I added Gaus- 
sian white noise (with signal-to-noise ratios (SNR) of 100, 10, or 2) to the hill images. 
The resulting brightness values were then clipped to the range [0, 1]. The results of 
running the z-only algorithm on these noisy images are shown in Figure 7-1 and Ta- 
ble 7-1. The figure shows that the algorithm performs well even for a signal-to-noise 
ratios of 2 and that the estimate degrades gracefully as the SNR is decreased. 

7.2 Camera Geometry Error 

The effect of errors in the camera geometry were investigated by adding error to the 
baseline distance, interior orientation parameters, or by rotating the images to move 
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Case 


Change 


Rel. Error 


Abs. Error 




Nominal (z-only) 


0.0175 


0.0247 


1 


SNR 100 


0.0203 


0.0275 


2 


SNR10 


0.0751 


0.1323 


3 


SNR 2 


0.1436 


0.1701 



Table 7-1: RMS Estimated surface error from noisy images. 
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Figure 7-1: Performance of z-only algorithm on noisy hill images 



the epipolar lines out of alignment. Each change was introduced independently. The 
baseline error images were created by changing both the baseline and nominal depth 
by the same percentage in order to keep the interior orientation parameters constant 
(10%, 20%, and 30% changes where introduced). The interior orientation error images 
were created by offsetting the origin of the camera coordinates in a random direction 
with magnitude of 1, 3, or 6 pixels. The epipolar error images were created by rotating 
the images so the total error introduced between the images was 0.5, 1, or 3 degrees. 
The results of running the z-only algorithm on these changed images are shown in 
Figures 7-2-7-4 and Table 7-2. The baseline errors mostly affect the absolute depth 
of surface estimate and have little effect on the relative error of the surface. On the 
other hand, the small change I introduced into the interior orientation parameters 
results in a profound change in the estimated surface (see Figure 7-3). Clearly it 
is important to have accurate knowledge of the camera parameters to obtain good 
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Figure 7-2: Performance of ^-only algorithm with hill baseline errors. 
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Figure 7-3: Performance of ^-only algorithm with hill interior orientation errors. 
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Case 



Change 



Rel. Error Abs. Error 





Nominal (z-only) 


0.0175 


0.0247 


1 


baseline error ( — 10%) 


0.0189 


10.0075 


2 


baseline error (—20%) 


0.0374 


99.9331 


3 


baseline error (—30%) 


0.0609 


199.8330 



1 


interior orientation error (1 pixel) 


0.0716 


0.0744 


2 


interior orientation error (3 pixels) 


0.3214 


0.3228 


3 


interior orientation error (6 pixels) 


0.9985 


1.0506 



1 


epipolar error (0.5°) 


0.0319 


0.0499 


2 


epipolar error (1°) 


0.0570 


0.0630 


3 


epipolar error (3°) 


0.0943 


0.2432 



Table 7-2: RMS Estimated surface error from geometry errors. 

estimates. 

Figure 7-4 shows that a moderate amount of error in the epipolar calibration can 
be tolerated with only a small effect on the estimated surface. This is good news 
since most planetary images will have to be re-projected into the aligned optical axes 
coordinate system. Figure 7-4 shows that the algorithm can tolerate some errors in 
this projection and still produce meaningful results. 



7.3 Reflectance map errors 

The effects of errors in the reflectance map were investigated by changing the light 
source positions or by scaling brightness values in the images to simulate an error 
in albedo calibration. I did not investigate the effects of using the wrong reflectance 
(such as Minnaert reflectance) since the estimation errors due to this effect cannot be 
generalized. The light source error images were created by changing the light source 
direction by 5, 15, or 30 degrees in a random direction. The albedo error images were 
created by scaling the images by 1/0.99, 1/0.95, or 1/0.90. The results of running the 
z-only algorithm on these changed images are shown in Figures 7-5-7-6 and Table 7-3 
Figure 7-5 shows that light source position errors of up to 30 degrees have relatively 
little effect on the estimated surface. While, typical photo-topography images have 
light source (sun) positions that are known to high precision, these results show that 
the z-only algorithm could be applied to images where the light sources are not as 
well known. 

On the other hand, Figure 7-6 shows that the algorithm can tolerate albedo cali- 
bration errors up to 10%. These results are misleading, however, since the algorithms 
fail to converge when the albedo error is greater than 10%. The algorithm may not 
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Figure 7-4: Performance of ^-only algorithm with hill epipolar errors. 
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Figure 7-5: Performance of z-only algorithm with hill light source errors. 
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Case 



Change 



Rel. Error Abs. Error 





Nominal (z-only) 


0.0175 


0.0247 


1 


light source error (5°) 


0.0551 


0.1912 


2 


light source error (15°) 


0.0178 


0.3219 


3 


light source error (30°) 


0.0922 


1.1150 



1 


albedo error (1%) 


0.0241 


0.0341 


2 


albedo error (5%) 


0.0681 


0.0686 


3 


albedo error (10%) 


0.1234 


0.1253 



Table 7-3: RMS Estimated surface error from reflectance errors. 

converge in those cases since the Lambertian reflectance map cannot exceed 1.0 and 
the error images produced contain many brightness values beyond that limit. The al- 
gorithm can probably tolerate albedo calibration errors that don't create normalized 
brightness values greater than 1.0. 

7.4 Summary 



The foregoing plots and tables indicate the z-only algorithm is fairly robust. In 
particular, they indicate that the algorithm can produce a reasonably good estimate 
even if the images are noisy, the camera geometry is not known perfectly, and the 
reflectance properties are in error. They also show that it is very important to have 
accurate internal orientation parameters, and an accurate baseline in order to estimate 
the depth correctly. 
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Figure 7-6: Performance of ^-only algorithm with hill albedo errors. 
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Chapter 8 
Algorithm Extensions 



In this chapter I discuss two ways of extending the algorithm to varying albedo 
surfaces. While all the previous algorithms have been restricted to surfaces with a 
constant (known) albedo, the algorithms discussed in this chapter work for surfaces 
that have markings or striations. The new algorithms still require that the geometric 
reflectance properties be constant and known for the whole surface. 

I also show in this chapter how more general camera geometries can be accommo- 
dated. Basically, the images are projected into a coordinate system that has aligned 
optical axes. 

8.1 Varying albedo algorithms 

As mentioned in Section 2.5, the simplification of constant albedo severely restricts 
the applicability of the algorithms that are developed in Chapter 4. In this section, I 
lift that restriction. The algorithms that result, do converge to a solution close to the 
actual surface, but the convergence is slower than for the constant albedo algorithms. 

The hrst thing to understand is whether we would expect a varying albedo algo- 
rithm to work. In other words, 'Does it seem reasonable that a unique value of albedo 
can be chosen for each point in the image'? Let's investigate that question. 

Each point in the images provides two constraints via the two extended photo- 
topography equations 

E {1) (x + ^,y) = p(x,y)RW(p,q), 

fb (8 - 1} 

EW(x - J —,y) = p(x,y)RW(p,q), 

for the two unknowns z and p. Since this is a situation with two equations and 
unknowns, a solution is at least conceivable. Investigating further we find that for 
a given p(x } y) there are at most two gradient directions (p(x } y) } q(x } y)) that can 
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satisfy each equation. 1 These gradient directions cannot be chosen arbitrarily since 
they must be consistent with the underlying height z, via the integrability constraint 
equations, 

fZx 

P = — , 

X J + Z (8.2) 

q y 



yz y + z 

A solution is obtained when a common, consistent gradient estimates exist. Con- 
versely, a given z defines the gradient components and the relationship between points 
in the images. A solution is obtained when p can be chosen to match the images E^' 
with the images estimated from the reflectance maps. 

From the foregoing discussion it should be clear that it is possible to create image 
sets that are inconsistent so that no solution exists. On the other hand it should 
also be clear that given consistent images a solution exists. Thus we find that the 
extended problem is well-posed. 

8.2 Minimizing departure from a constant albedo. 

The hrst algorithm is applicable to slowly varying albedo surfaces. It is based on the 
z-only cost function but includes a penalty term of the form p(p — p) 2 where p is 
the average albedo over the whole image. This additional term penalizes departure 
from a constant albedo and can be used to estimate the calibration factor Xp on the 
reflectance. The cost function including this term is 



mm 



' mJ= lll{ { E(1){X + 2? y) ~ P{X ' y)R(1){P ' q) ) 

+ (E^(x-^,y)-p(x,y)^{p,q) > j 

+ A [zl + 2z 2 xy + 4j + p(p - p) 2 } dx dy (8.3) 

where 



J v pdx dy 



(8.4) 



Iv dx d v 

and T> is the whole image. The smoothness penalty term is used to guide convergence 
and is not required to guarantee a unique solution. Usually, A is slowly reduced to 
zero as the solution is reached to avoid biasing the solution. The algorithm based 
on this cost function converges, albeit slowly, and works best for surfaces that have 
nearly constant albedo. 



L For surfaces with Lambertian Refiecta 



8.3. MINIMIZING LOCAL ALBEDO CHANGE. 
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8.3 Minimizing local albedo change. 

The second algorithm is applicable to surfaces with piecewise constant or piecewise 
linear albedo. In this case, the cost function includes a penalty term of the form 
p(p — p) 2 where p is the average in some local neighborhood J\f. This term is similar 
to a discrete approximation to the Laplacian and penalizes departure from a local av- 
erage. For a piecewise constant albedo, this term will be zero except on the boundary 
between the constant areas of albedo. The cost function including this term is 



min J 



where 



+ (e^(x - ^,y) - p(x,y)R^(p,q)\ 
+ A [z 2 xx + 2z 2 xy + 4j + p(p - p) 2 } dx dy 

_ $ M pdx dy 



(8.5) 



(8.6) 



J M dx dy ' 
Following the implementation of the z-only algorithm, the cost function can be 



discretized to become, 
min J 



1 



2NMe 2 



x,yE.V 



E (^ 



(x,y) 



p{x,y 



)R {1) (p, 



+ {FW(x,y)-p(x,y)RW(p,q) 

+ X \ Z lx + 2z ly + Z yy] + »(P ~ P) 2 } 

where p is the average in a 3-by-3 neighborhood, 



(8.7) 
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The performance of this algorithm on the varying albedo test images is shown in 
Figures 8-2-8-4. 



112 CHAPTER 8. ALGORITHM EXTENSIONS 

8.4 Varying albedo test results. 

To test the algorithm based on this cost function, two of the test images were modified 
to include either a p = 0.7 albedo strip from the lower left to the upper right (Figure 8- 
1) or an p = 0.7 albedo variation in a layer similar to what would be expected from 
sedimentation processes (Figure 8-3). 

8.4.1 Crater with dark stripe 

The hrst set of test images has the p = 0.7 albedo stripe imposed on the easy crater 
test case. The light source geometry, true surface, and true albedo are shown in 
Figure 8-1. The albedo variation is "painted" on the crater surface. 

The Figure 8-2 shows the result of applying the varying albedo algorithm to this 
test case. The figure shows the cost function history and optimization parameters as a 
function of the number of function evaluations. Also shown are the estimated surface 
and albedo image at various stages during the convergence. While the algorithm 
does converge to a surface and albedo map close to the true values, the convergence 
is much slower than for the constant albedo algorithm and the final surface estimate 
has more error. Of course, the constant albedo algorithm would not perform any 
better if presented with the varying albedo images. 

8.4.2 Hill with sedimentation 

The other set of varying albedo test images is based on the hill test case with the 
albedo set to p = 0.7 for surface heights in a certain range (see Figure 8-3). The light 
source positions and geometry are the same as for the hill test case. 

Figure 8-4 shows the result of applying the varying albedo algorithm to this test 
case. The figure shows the cost function history and optimization parameters as a 
function of the number of function evaluations. Also shown is the estimated surface 
and albedo image at various stages during the convergence. While the algorithm 
does converge to a surface and albedo map close to the true values, the convergence 
is slower than for the constant albedo algorithm and the final surface estimate has 
more error. 

8.4.3 Summary 

The varying albedo algorithm can successfully estimate both the surface depth and 
the albedo variation but requires at least 2 times as many function evaluations than 
the constant albedo z-only algorithm. In fact, the performance figures in this section 
show surfaces that have not converged even after 2800 function evaluations. However, 
these algorithms must be used for surfaces that cannot be estimated correctly with 
the z-only algorithm because of varying albedo. 
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Figure 8-1: Crater on flat plane with varying albedo. Shown is the camera geometry as 
projected into the xz- and y^-planes, the true surface as a mesh plot, albedo image, and 
reflectance function contours for the two light source positions. 
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Figure 8-2: Performance of the varying albedo algorithm on the varying albedo crater 
images. 
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Figure 8-3: Test images of hill with varying albedo. Shown is the camera geometry as 
projected into the xz- and y^-planes, the true surface as a mesh plot, albedo image, and 
reflectance function contours for the two light source positions. 
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Figure 8-4: Performance of the varying albedo algorithm on the varying albedo hill images. 
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Figure 8-5: Camera reprojection geometry. 

The slow convergence is a direct result ol increasing the number ol degrees ol 
freedom available to the algorithm. It has the choice of meeting the photo-topography 
constraints by changing the depth or by just changing the albedo. Initially it's easier 
to change the albedo since it shows up directly in the constraint equations. However, 
using albedo changes alone is not sufficient to minimize the equations, so after a while 
the surface height begins the change. 

The convergence would probably be faster if the stereo part of the algorithm was 
more powerful. The disparity of the albedo edges should provide a strong constraint 
on the height. The fact that they don't is due to the weakness of the stereo part of the 
algorithm. Hopefully ongoing research will turn up a more powerful stereo algorithm. 



8.5 General Camera Geometry 

The algorithms presented thus far are based on the simplification that the camera 
optical axes are aligned. This is never true for images taken from moving vehicles 
such as inner-planetary probes, satellites, or aircraft. The simplification of aligned 
optical axes is used since it greatly simplified the equations. The simplification can 
be lifted by re-projecting the images from their true camera coordinates into the 
coordinates of a virtual camera with its principal point in the same position but 
with the aligned orientation (see Figure 8-5). I refer to the true camera coordinates 
as the primed coordinate system (i.e., (x',y',z') T and the aligned coordinates using 
un-primed notation, consistent with the equations presented so far. 

Points in a camera image map to rays in space. The mapping between coordinate 
systems preserves the orientation of these rays in the global coordinate system. Fet 
r' = (x', y', f') T be the coordinates of a point in the camera's image. This maps to 
the ray in the direction r'. Fet T be the rotation matrix relating the two coordinate 
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systems, that is, 

r' = Tr, (8.9) 

for any vector r (see Figure 8-5). The ray in the aligned coordinate system is then 
along the direction Tr 1 . This ray can be re-projected into the virtual image plane 
with focal length / by normalizing, 

- - /(TV » (8.10) 



(T T r') • z 

where z = (0,0, 1) T is the unit vector along the z axis in the aligned coordinate 
system. Similarly, 

r'=7^l (8.1D 

(Tr) • z' y J 

relates points in the un-primed image to points in the primed image. 

We can use these equations to create the image that would have been seen by a 
camera with aligned optical axes. Given E'(r') in the true camera coordinate system, 
the re-projected image is 

£(r) = e {$$*) (8 - 12 » 

since T T is the rotation matrix from the aligned coordinates to the true coordinates. 
For r = (x, y, f) T this equation written in (x, y) notation is 

cv x rpiMiX + t 12 y + W t 21 x + t 22 y + t 23 f 

E{x,y) = E{- — — r-7,T ~7f 777, (8-13) 

hix + t 32 y + t 33 f t 31 x + l 32 y + t 33 j 

where the elements of the 3-by-3 matrix T are tij. For discrete images, the value of 
the above expression can be computed using some type of interpolation (say bilinear 
interpolation). 

Note that this re-projection can be easily incorporated into the z-only, disparity- 
based, and varying albedo algorithms by defining F^' to interpolate in E' instead of 
E. 



Chapter 9 

Summary 



This thesis has presented a methodology for combining or fusing multiple vision algo- 
rithms. Four different cost functions (and their associated algorithms) were presented 
that illustrate the methodology. The basic methodology is to combine the constraint 
equations of the problem to form a single cost function in the spirit of variational 
calculus. 

The performance of the four algorithms was evaluated using four synthetic noise- 
free test images of varying difficulty The most closely-coupled algorithm, the z-only 
algorithm, had the best performance. The z-only algorithm was able to correctly 
estimate the synthetic surface in about 200 function evaluations for three of the four 
cases. For the remaining case (the hard crater images), the algorithm got stuck in a 
local minimum, as did all the algorithms tried. This case has lighting geometry that 
results in ambiguous shading information. 

It was shown that the z-only algorithm has much better performance using the 
two photo-topography images than a simple shape- from-shading algorithm which uses 
only one image. This performance increase validates the fusion approach to obtaining 
better performing vision algorithms. 

The z-only algorithm was also shown to be robust; able to accurately estimate the 
synthetic surface in the presence of several types of errors. The performance of the 
algorithm based on images that contained noise, geometry error, or reflectance errors 
was shown. In most cases, the algorithm was able to form a good estimate. 

9.1 Mars Images 

The robustness and performance of the algorithm on synthetic images, builds confi- 
dence that the algorithm can perform similarly on real images. One such set of real 
images is shown in Figure 9- 1. 1 The images are Viking stereo images of Mars. The 
images as received (and shown) are processed versions of the original Viking images. 



1 We are grateful to Mike Caplinger for providing these images. 
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They have been reprojected into an aligned coordinate system and hltered to remove 
biases due to any large scale albedo variations. The images also have very similar 
lighting conditions (they were probably taken close together in time). This image 
pair thus doesn't represent the best possible situation for the algorithms developed 
in this thesis. Nevertheless, the z-only algorithm performs well. 

The estimated surface from the z-only algorithm is shown in Figure 9-2 and 9- 
3. The figures show that the z-only algorithm produces a reasonable estimate 
of the Martian surface. It should be noted that the camera geometry and light 
source positions were given, but the reflectance map for the surface of Mars was not. 
The results shown are based on the assumption of a Lambertian reflectance map. 
The nominal slope of the estimated surface (see surface plots) is an artifact of the 
camera baseline orientation with respect to the planet's surface and can be removed 
by shifting the camera positions along their line of sight. In fact, the estimate shown 
is based on cameras that are shifted slightly so that the light source positions could 
be represented using (p S} q s ) gradient components. 2 

While it is difficult to see in the halftone reproduction of the estimated images, 
the images contain a slight ghosting effect. 3 The ghosting is like a double image and 
is depth dependent. By running the algorithm with slightly shifted camera principal 
points, different parts of the estimated images can be brought into registration. The 
ghosting is probably caused by using Lambertian reflectance instead of the true radi- 
ance function for the Martian surface. However, even though the algorithm used an 
inaccurate reflectance function, the surface estimate does contain many of the small 
and large scale features found in the Viking images. Examples include the cliff in the 
upper right and the valley just below, as well as the craters near the bottom of the 
images and the ridges in the center. An even better estimate of the Martian surface 
could have been obtained if the true reflectance function were used with the z-only 
algorithm. 

9.2 Future research 

There are several directions this research can go in the future: 

• Apply the fusion methodology outlined in this thesis to other vision problems. 

• Apply the z-only algorithm to more real images to further judge its performance. 
Particularly apply the algorithm to image sets with differing lighting conditions. 

• Determine the performance of the z-only algorithm when more realistic re- 
flectance maps are used with the algorithm. Recall that while I only tested the 



2 A requirement for the particular implementation I have used. With slight modifications, the 
algorithms could be used with any other light source position parameterization. 
3 The extent of the ghosting is less than 2 pixels for a 257-by-257 image. 



9.2. FUTURE RESEARCH 



121 



Camera Geometry 



500 



-500 



a- -looo 



-1500 



-2000 



c^ 



cM^ 




500 



Reflectance Function Contours 



cr 





True images 



Figure 9-1: Stereo images of Mars taken by the Viking probe. 
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Figure 9-2: Results of running the z-only algorithm on the Viking Mars images. 
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Figure 9-3: High resolution mesh and surface plot of estimated Mars surface. The surface 
plot was created by coloring the surface using the left estimated image. 
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algorithms on surfaces with Lambertian reflectance, other reflectance lunctions 
can easily be used with the algorithm. 

• Develop an algorithm with better stereo integration that can correctly estimate 
the surface for the Hard Crater test case. To do this may require developing a 
new innovative stereo algorithm. 

Of these possible directions, the last one presents the most challenges and offers the 
most rewards. 

9.3 Conclusions 

This thesis has shown that the variational approach to fusion problems can produce 
robust, well performing vision algorithms. It has also shown that fused algorithms 
can have significant performance gains over non-fused alternatives. 

The four example algorithms, each based on a different variable representation, 
showed that choosing the right variable representation is important to achieving good 
estimates. 

Finally in conjunction with this research, a new algorithm, the z-only algorithm, 
was developed to solve the photo-topography problem. This is the first well perform- 
ing algorithm to solve that problem. 



Appendix A 

Gradient derivation for ^-only 
algorithm 



This appendix presents the gradient derivation for the z-only algorithm for use with 
the conjugate gradient optimization technique. 

A.l Cost function. 

In the main text, the discrete approximation to the cost function for the z-only 
algorithm is given as 

min „ 

2MNe c _ 



J = ^^ E { i F{1) ^ y) ~ R{1) (P' ?)) 2 + i F{2) ^ y) ~ R{2) (P' *)) 

+ \[z 2 xx + 2z 2 xy + z 2 yy \} (A.r 



where T> is the discrete domain of the underlying variables in the global coordinate 
system, p and q are computed using 

fZx 

p 



XZt , + Z (A.2) 

yz y + z 

M and N are the row and column dimensions of the discrete domain and e is the 
grid spacing (assumed to the same in both the x and y directions). The F^'(x,y) are 
interpolated from the input images E^'(x,y) via linear interpolation, 

F (i \x,y) = E (i \x,y) + (x ± ^ - x) [E (i \x + l,y)- E (i \x,y)\ (A.3) 
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where 

r i 

x =floor(x± l~). (A.4) 

and the floor(x) function returns the greatest integer that is smaller than x. These 
equations assume that x is sampled on unity-spaced grid. 

A. 2 Bicubic extrapolation 

Recall that the matched-grid implementation is used for the z-only algorithm. In 
this implementation, z is the same size as the image arrays E^' and the boundary is 
extended if necessary based on bicubic interpolation. Let zn be a value in the hrst 
column of z (i.e., on the left boundary of z). Then the extrapolated value z j0 , is 
computed using 

z i0 = oZn — OZi2 + Zi 3 . (A. 5) 

Similar equations work for the other boundaries of z. In this way, a M-hy-N matrix 
can be extrapolated to form a (M + 2)-by-(iV + 2) matrix. 

A. 3 Gradient of terms based on convolution. 

Theorem 1 Let y be computed via a 2-D convolution from the field z and the filter 
h, y = z * h. Let f(x) : 5ft — > 9? ; be any scalar point function. Then the derivative of 
the scalar cost function J = Y^iY^j f(yij) with respect to z\ m is 



dJ df(x) 



OZ^ rn CIX 

where 



*h (A.6) 



h l3 = h-i-j. (A. 7) 
Proof Using the definition of the 2-D convolution, 

Vij =^2^2 Z kshi-k,j-s, (A. 8) 

k s 



the cost function can be expanded to obtain, 



i j \ k s 



i — k,j — s 



(A.9) 
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The sums are taken over the entire range of indices into the matrices. Taking the 
derivative and using the chain rule we find, 



d J Tr-.sr-.dfd 



dzi 



* 3 



dx dzi Ti 



k s 



^^dx^^ dzi~ ,j ~ s 

* 3 

* 3 



df_ 
dx 



* h 



(A.10) 
(A.11) 

(A.12) 
(A.13) 
(A.14) 



The result holds for computational molecules as well since operating with a com- 
putational molecule m is equivalent to convolving with rh. 

Theorem 2 Theorem 1 also holds for scalar functions of multiple arguments. In that 
case, 

d.J df 7 df 



oz ax\ ox 2 



(A.15) 



Proof Suppose y and w are defined via 2-D convolution, y = z * h, w = z * g, 
and f(xi 7 x 2 ) : [9£ X 9£] — > 9£, then the partial derivative of the cost function J = 
T,iT,j f(yij,Wij) is 

"^~ = ^~~ * h + ^~ *9- (A.16) 

OZ OX\ OX 2 

The proof follows the same lines as the proof for Theorem 1 I. 



A. 4 Cost function derivative. 

Using the chain rule, the derivative of the cost function (Equation A.l) is 



dJ 



dz 



km 



£W(FW-W 

x y I \ 

+ {F (2)_ R (2)J 



(JZ km 

ap(i) 









dz km p dz k 



oz kn 



+ A 



'J% xx UZ X y 'JZyy 

%t.t. ~ \ ^Z X y — -\- Z t 



dz kri 



dz km dz k 



(A17) 
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Let's look at each term in turn. The derivatives involving the interpolated images 
can be computed using the definition of F^\ 

dF E x %y)p- + UP, ~ ^-\ [E(x + l,y) - E(x,y)] 



OZ km OZj~ m y LZ OZkr, 

+ (x±^--x) [E x (x + l,y)- E x (x, y)] -^- (A.f 8) 



11 

2z"< 



T^E x ((x),y). (A.f 9) 



That last equality is possible since dx/dz^m = for one-sided derivatives. 

The derivatives of the gradient components p and q can also be computed using 
their definition. They are, 



dz x / „ \ „\ f~(„ dz x i dz 



dp f&(*** + *)~f **(*& + 



km ^^km ^^km 



£'~t 



anc 



(JZfc m 


(xz x + z) 2 




f„ dz x r„ dz 

° z km c ''km 




(xz x + z) 2 


dq 


- 1 dz km J ^y dz km 



dz km (yz y + z) 



(A.20) 
(A.21) 

(A.22) 



As discussed in the main text, the partial derivatives of z (z X} z y} z XX} and so on), 
are computed via 2-D computational molecules. Call these molecules, h X} h y} h xx , 
and so on, in the obvious way. Using the results of Theorem f , the derivative of the 
cost function can then be written, 

dJ r -(F^ - R^)E^ + (F^ - B?))EW] fb 



OZ^'m ZZ 

f z v 



+ [(FW - R^ + (i*> - R^} ( -^- 



+ [( F (D _ rW)rW + (F^ - rW)rW] 



yZy + Z) 



(FU - R^)RV + (F^ - B?))B£)] . f " - ) o K 

K J (xz x + z) 2 ) 



T A LZ^a; O /l^j; + ZZ X y O fl X y -\~ Zyy O tlyy\ ■ ( A . Z 6 J 

These equations are not the whole story since they don't work on the boundary. I 
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will show how to take into account the effect of the bicubic extrapolation on the left 
boundary; the other boundaries are dealt with in a similar way. Suppose h has either 
2 or 3 columns, and let z be the extrapolated version of z (for simplicity assume z 
is square of size JV-by-iV). In this case the subscripts on z run from to N — 1. 
Theorem 1 states that the derivative of J = Y, J2 f(z * h) is 

* = &•*• 

Call this matrix G. Since z is computed from z via the bicubic extrapolation described 
in Section A. 2, we can compute dJ/dz by taking into account the dependence of z^ 
on zn, Zi 2 , and Z&. Let G be the matrix dJ/dz, then columns 1-3 of G can be 
computed using 

Gi\ = Gi\ + 3G 8 o, (A. 25) 

G i2 = G i2 - 3G i0 , (A.26) 

G i3 = G i3 + G i0 . (A.27) 

Similar expressions hold for the other boundaries. 
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Appendix B 
M-file Listings 



This appendix contains the MATLAB 1 M-file source hies for the DFSS algorithms. 
The source hies are presented for the six algorithms in this thesis: zpq } z-only, dual 
z, disparity, sfs, and varying albedo. A key to the hies is given in Table B-l. 

B.l Cost function and gradient routines 
B.l.l hbdfss_cost3c.m 

l function [J ,FR1 , FR2]=hbdf ss_cost (v,zsize , levels , params , El , E2 , lambda, mu) 



12 
13 
14 
15 
16 
17 
18 



HBDFSS_C0ST3C Cost function for depth from shading and stereo 
problem using hierarchical basis functions. Uses true 
perspective projection. 

J=HBDFSS_C0ST3C (V,ZSIZE, LEVELS, PARAHS, El ,E2, LAMBDA, HU) where 

V = [z,p,q] are the optimization variables, ZSIZE is 

SIZE(z), LEVELS is the number of h-basis levels. El and E2 are 

the input images, LAMBDA is the scalar weighting factor on 

departure from smoothness and MU is the scalar weighting factor on 

integrability . The image parameters 

PARAMS = [f ,b,zO,psl ,qsl ,ps2 ,ps2 ,vxl , vyl ,vx2,vy2] . 

Clay M. Thompson 5-18-92 
Revised to use correct p,q calculation. 
Revised to use p,q,z of the same size as E. 
Revised to output estimated image. 
Revised to support multi-grid scheme. 



20 '/, Camera constants 

21 f = params (1) ; 

22 b = params (2) ; 

23 zO = params(3); 

24 gamma = f *b/2 ; 

25 delta = params(12); 

26 area = prod(zsize)*delta*delta; 

27 

28 '/, Light Source positions 

29 psl = params(4); qsl = params(5); 

30 ps2 = params(6); qs2 = params(7); 

31 

32 '/, Camera coordinate calibration 

33 vxl = params(8); vyl = params(9); 

34 vx2 = params(lO); vy2 = params(ll); 



1 Matlab is a product of The MathWorks, Inc., 24 Prime Park Way, Cochituate Place, Natick 
MA 01760, (508)653-1415. 
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Algorithm cost functions and gradient routines 



hbdfss_cost3c.m 

hbdfss_grad3c.m 

hbdfss_cost2c.m 

hbdfss_grad2c.m 

hbdfss_cost4.m 

hbdfss_grad4.m 

hbdfss_cost7.m 

hbdfss_grad7.m 

hbdfss_cost8b.m 

hb dfss_grad8b .m 

hbsfs_cost.m 

hbsfs_grad.m 



Cost function for 
Gradient function 
Cost function for 
Gradient function 
Cost function for 
Gradient function 
Cost function for 
Gradient function 
Cost function for 
Gradient function 
Cost function for 
Gradient function 



zpq algorithm. 

for zpq algorithm. 
z-only algorithm. 

for z-only algorithm, 
dual- z algorithm. 

for dual- z algorithm, 
disparity algorithm. 

for disparity algorithm, 
varying albedo algorithm. 

for varying albedo algorithm, 
shape-from-shading algorithm. 

for shape-from-shading algorithm. 



Support routines 



rmap Lambertian reflectance function. 

rmapp Derivative of reflectance function with respect to p. 

rmapq Derivative of reflectance function with respect to q. 

conjgrad.m Conjugate gradient optimization. 

lsearch.m Line search function for conjugate gradient optimization. 

filter2d.m 2-D computational molecule filtering. 

cfilter2d.m 2-D computational molecule filtering with bicubic interpolation. 

hbasis.m Main level hierarchical basis conversion. 

hb.m Hierarchical basis interpolation. 

hbt.m Adjoint hierarchical basis interpolation. 

interpx.m Linear interpolation in the x direction. 

domain2d.m 2-D plaid domain generation. 

icubic.m 1-D cubic interpolation. 

dcubicx.m Derivative of 1-D cubic interpolation with respect to x. 

dcubicz.m Derivative of 1-D cubic interpolation with respect to z. 



Table B-l: M-file Descriptions. 
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35 vx = (vxl+vx2)/2; vy = (vyl+vy2)/2; 

36 

37 '/, Extract z,p,q and transform into nodal basis. 

38 mz = zsize(l); nz = zsize(2); 



,l:nz); 

, nz+[l :nz] ) ; 

,2*nz + [1 :nz] ) ; 
) = hbasis(z , levels) ; 
) = hbasis(p, levels) ; 
) = hbasis(q, levels) ; 



39 Z = v( 

40 p = v( 

41 q = v( 

42 z( 

43 p( 

44 q( 
45 

46 '/, Spacial coordinates in image. 

47 [mz,nz] = size(z); 

48 [x,y] = domain2d([0:mz-l]*delta+vx+l , [0 :nz-l]*delta+vy+l) ; 

49 

50 '/, Stencils 

51 hxz = [-1 l]/(2*delta) ; '/, X-derivative for zx 

52 hyz = [1 ;0;-l]/(2*delta) ; '/, Y-derivative for zy 

53 hx = [-1 1;-1 l]/(2*delta) ; '/. X-derivative 

54 hy = [1 1;-1 -l]/(2*delta) ; '/. Y-derivative 

55 

56 '/, Compute reflectance map values. 

57 Rl = rmap(p,q,psl ,qsl) ; 

58 R2 = rmap(p,q,ps2 ,qs2) ; 

59 

60 '/, Compute numerical derivatives of p,q, and z. 

61 px = f ilter2d(p,hx, 'resize ') ; 

62 py = f ilter2d(p,hy , 'resize ') ; 

63 qx = f ilter2d(q,hx, 'resize ') ; 

64 qy = f ilter2d(q,hy , 'resize ') ; 

65 zx = cf ilter2d(z ,hxz) ; 

66 zy = cf ilter2d(z ,hyz) ; 

67 

68 '/, Compute estimates of p,q using numerical derivatives of z. 

69 pe = f *zx. /(x . *zx+z) ; 

70 qe = f *zy . /(y . *zy+z) ; 
71 

72 '/, Compute disparity. 

73 d = (f*b/2)./z; 
74 

75 '/, Determine stereo mapped images F1,F2. Set error to zero where F is NaN 

76 Fl = interpx(El, (x+d-vxl-l)/delta+l) ; 

77 out = isnan(Fl); if any(out ( : )) , Fl(out) = Rl(out); end, 

78 F2 = interpx(E2, (x-d-vx2-l)/delta+l) ; 

79 out = isnan(F2); if any(out ( : )) , F2(out) = R2(out); end, 

80 

81 terml = (0 . 5/area)*sum(sum((Fl-Rl) . ~2 + (F2-R2) . ~2) ) ; 

82 term2 = (0 . 5/area)*lambda*sum(sum(px . ~2 + py.~2 + qx.~2 + qy.~2)); 

83 term3 = (0 . 5/area)*mu*sum(sum((pe-p) . ~2 + (qe-q).~2)); 

84 '/.disp(sprintf ( 'Terms: '/.12.5f '/.12.5f '/.12 .5f ', terml ,term2 ,term3) ) ; 

85 J = [terml+term2+term3 , terml ,term2 ,term3] ; 

86 

87 if nargout>l, 

88 FR1 = [El;Fl;Rl;abs(Fl-RD] ; 

89 FR2 = [E2;F2;R2;abs(F2-R2)] ; 

90 end 

91 



B.1.2 hbdfss_grad3c.m 



l function [J ,FR1 , FR2]=hbdf ss_cost (v,zsize , levels ,params , El , E2 , lambda, mu) 
HBDFSS_C0ST3C Cost function for depth from shading and stereo 
problem using hierarchical basis functions. Uses true 
perspective projection. 

J=HBDFSS_C0ST3C(V,ZSIZE, LEVELS ,PARAMS, El ,E2 , LAMBDA, MU) where 

V = [z,p,q] are the optimization variables, ZSIZE is 

SIZE(z), LEVELS is the number of h-basis levels. El and E2 are 

the input images, LAMBDA is the scalar weighting factor on 

departure from smoothness and MU is the scalar weighting factor on 

integrability . The image parameters 

PARAMS = [f ,b,z0,psl ,qsl ,ps2 ,ps2 ,vxl , vyl ,vx2,vy2] . 

Clay M. Thompson 5-18-92 
Revised to use correct p,q calculation. 
Revised to use p,q,z of the same size as E. 
Revised to output estimated image. 
Revised to support multi-grid scheme. 



12 
13 
14 
15 
16 
17 
18 
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20 '/, Camera constants 

21 f = params(l) ; 

22 b = params (2) ; 

23 zO = params(3); 

24 gamma = f *b/2 ; 

25 delta = params(12); 

26 area = prod(zsize)*delta*delta; 

27 

28 '/, Light Source positions 

29 psl = params(4); qsl = params(5); 

30 ps2 = params(6); qs2 = params(7); 

31 

32 '/, Camera coordinate calibration 

33 vxl = params(8); vyl = params(9); 

34 vx2 = params(lO); vy2 = params(ll); 

35 vx = (vxl+vx2)/2; vy = (vyl+vy2)/2; 

36 

37 '/, Extract z,p,q and transform into nodal basis. 

38 mz = zsize(l); nz = zsize(2); 

39 z = v( : ,1 :nz) ; 

40 p = v( : ,nz+[l :nz] ) ; 

41 q = v(:,2*nz + [l:nz]); 

42 z(:) = hbasis(z , levels) ; 

43 p(:) = hbasis(p, levels) ; 
) = hbasis(q, levels) ; 



44 q( 

45 



46 '/, Spacial coordinates in image. 

47 [mz,nz] = size(z); 

48 [x,y] = domain2d([0:mz-l]*delta+vx+l , [0 :nz-l]*delta+vy+l) ; 

49 

50 '/. Stencils 

51 hxz = [-1 l]/(2*delta) ; '/. X-derivative for zx 

52 hyz = [1 ;0;-l]/(2*delta) ; '/. Y-derivative for zy 

53 hx = [-1 1;-1 l]/(2*delta) ; '/. X-derivative 

54 hy = [1 1;-1 -l]/(2*delta) ; '/. Y-derivative 

55 

56 '/, Compute reflectance map values. 

57 Rl = rmap(p,q,psl ,qsl) ; 

58 R2 = rmap(p,q,ps2 ,qs2) ; 

59 

60 '/, Compute numerical derivatives of p,q, and z. 

61 px = f ilter2d(p,hx, 'resize ') ; 

62 py = f ilter2d(p,hy , 'resize ') ; 

63 qx = f ilter2d(q,hx, 'resize ') ; 

64 qy = f ilter2d(q,hy , 'resize ') ; 

65 zx = cf ilter2d(z ,hxz) ; 

66 zy = cf ilter2d(z ,hyz) ; 

67 

68 '/, Compute estimates of p,q using numerical derivatives of z. 

69 pe = f *zx. /(x.*zx+z) ; 

70 qe = f *zy . /(y . *zy+z) ; 
71 

72 '/, Compute disparity. 

73 d = (f*b/2)./z; 
74 

75 '/, Determine stereo mapped images F1,F2. Set error to zero where F is NaN 

76 Fl = interpx(El, (x+d-vxl-l)/delta+l) ; 

77 out = isnan(Fl); if any(out ( : )) , Fl(out) = Rl(out); end, 

78 F2 = interpx(E2, (x-d-vx2-l)/delta+l) ; 

79 out = isnan(F2); if any(out ( : )) , F2(out) = R2(out); end, 

80 

81 terml = (0 . 5/area)*sum(sum((Fl-Rl) . ~2 + (F2-R2) . ~2) ) ; 

82 term2 = (0 . 5/area)*lambda*sum(sum(px . ~2 + py.~2 + qx.~2 + qy.~2)); 

83 term3 = (0 . 5/area)*mu*sum(sum((pe-p) . ~2 + (qe-q).~2)); 

84 '/.disp(sprintf ( 'Terms: '/.12.5f '/.12.5f '/.12 .5f ', terml ,term2 ,term3) ) ; 

85 J = [terml+term2+term3 , terml ,term2 ,term3] ; 

86 

87 if nargout>l, 

88 FR1 = [El;Fl;Rl;abs(Fl-RD] ; 

89 FR2 = [E2;F2;R2;abs(F2-R2)] ; 

90 end 

91 

B.1.3 hbdfss_cost2c.m 

1 function [J ,FR1 , FR2]=hbdf ss_cost (z , n , params , El , E2 , lambda) 

2 */.HBDFSS_C0ST2C Cost function for depth from shading and stereo 

3 '/, problem with heirarchical basis functions. Uses true 

4 '/, perspective projection. 
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12 
13 
14 
15 
16 
17 
18 
19 
20 
21 



J=HBDFSS_C0ST2C(Z,N, PARAMS, El, E2 , LAMBDA) where Z is the surface 
height, N is the number of H-basis levels.. El and E2 are the input 
images, LAMBDA is the scalar weighting factor on departure from 
smoothness. The image parameters 

PARAMS = [f ,b,z0,psl ,qsl ,ps2 ,qs2 , vxl , vyl , vx2 , vy2 , delta] . 

J= [Jtot , Jimage , Jsmooth] . 

Clay M. Thompson 2-24-91 
Revised 4-3-91 by cmt 
Revised 5-10-91 by cmt 
Revised 5-19-91 by cmt 

Revised 4-30-92 by cmt to support multi-grid scheme. 
Revised to output estimated image. 



Camera constants 

22 f = params (1) ; 

23 b = params (2) ; 

24 zO = params(3); 

25 gamma = f*b/2/z0; 

26 delta = params(12); 

27 area = prod(size(z) )*delta*delta; 

28 

29 '/, Light Source positions 

30 psl = params(4); qsl = params(5); 

31 ps2 = params(6); qs2 = params(7); 

32 

33 '/, Camera coordinate calibration 

34 vxl = params(8); vyl = params(9); 

35 vx2 = params(lO); vy2 = params(ll); 

36 vx = (vxl+vx2)/2; vy = (vyl+vy2)/2; 

37 

38 '/, Spacial coordinates in image. 

39 [mz,nz] = size(z); 

40 [x,y] = domain2d([0:mz-l]*delta+vx+l , [0 :nz-l]*delta+vy+l) ; 

41 

42 '/, Stencils 

43 hx = [-1 l]/(2*delta) ; '/. X-derivative 

44 hy = [l;0;-l]/(2*delta) ; '/. Y-derivative 

45 hxx = [1-2 1; 2-42; 1 -2 l]/(4*delta*delta) ; 

46 '/.hxy = [-1 1;0 0;1 -1] /(4*delta*delta) ; 

47 hxy = [-1 1;1 -l]/(delta*delta) ; 

48 hyy =[12 l;-2 -4 -2 ; 1 2 1] /(4*delta*delta) ; 

49 del2 =[14 1;4 -20 4;1 4 1] /(6*delta*delta) ; 

50 

51 '/, Compute p,q using numerical derivatives of z. 

52 z = hbasis(z,n); 

53 zx = cf ilter2d(z ,hx) ; 

54 zy = cf ilter2d(z ,hy) ; 

55 p = f *zx . /(x. *zx+z) ; 

56 q = f *zy ./(y .*zy+z) ; 

57 

58 '/, Compute reflectance map values. 

59 Rl = rmap(p,q,psl ,qsl) ; 

60 R2 = rmap(p,q,ps2 ,qs2) ; 

61 

62 '/, Compute numerical derivatives of z. 

63 zxx = cf ilter2d(z ,hxx) ; 

64 zxy = cf ilter2d(z ,hxy ) ; 

65 zyy = cf ilter2d(z ,hyy ) ; 

66 

67 '/, Compute disparity. 

68 d = (f*b/2)./z; 

69 

70 '/, Determine stereo mapped images F1,F2. Set error to zero where F is NaN 

71 Fl = interpx(El, (x+d-vxl-l)/delta+l) ; 

72 out = isnan(Fl); if any(out(:)), Fl(out) = Rl(out); end, 

73 F2 = interpx(E2, (x-d-vx2-l)/delta+l) ; 

74 out = isnan(F2); if any(out ( : )) , F2(out) = R2(out); end, 
75 

76 terml = (0 . 5/area)*sum(sum((Fl-Rl) . ~2 + (F2-R2) . ~2) ) ; 

77 term2 = (0 . 5/area)*lambda*(sum(sum(zxx . ~2) ) + 2*sum(sum(zxy . ~2) ) + sum(sum(zyy . "2) ) ) ; 

78 */,term2 = (0 .5/area)*lambda*sum(sum(delsqz . "2) ) ; 

79 '/.disp(sprintf ( 'Terms: '/.12.5f '/.12 . 5f ', terml ,term2) ) ; 

80 J = [terml+term2 , terml ,term2] ; 

81 

82 if nargout>l, 

83 FR1 = [El;Fl;Rl;abs(Fl-RD] ; 
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84 FR2 = [E2;F2;R2;abs(F2-R2)] ; 

85 end 

B.1.4 hbdfss_grad2c.m 



6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 



function G=hbdf ss_grad(z ,n,params ,E1 ,E2 , lambda) 
HBDFSS_GRAD2C Gradient function for depth from shading and stereo 
problem with heirarchical basis functions. Uses true perspective 
projection . 

G=HBDFSS_GRAD2C(Z,N, PARAMS, El, E2 , LAMBDA) where Z is the height 
map, N is the number of H-basis levels. El and E2 are the input 
images, LAMBDA is the scalar weighting factor on departure from 
smoothness. The image parameters PARAMS = [f ,b ,z0 ,psl ,qsl ,ps2 ,ps2] ; 

Clay M. Thompson 2-24-91 

Revised 4-3-91 by cmt . 

Revised 5-10-91 by cmt. 

Revised 5-19-91 by cmt. 

Revised 4-30-92 by cmt for multi-grid scheme. 

'/. 



Camera constants 

18 f = params (1) ; 

19 b = params (2) ; 

20 zO = params(3); 

21 gamma = f*b/2/z0; 

22 delta = params(12); 

23 area = prod(size(z) )*delta*delta; 

24 

25 '/, Light Source positions 

26 psl = params(4); qsl = params(5); 

27 ps2 = params(6); qs2 = params(7); 

28 

29 '/, Camera coordinate calibration 

30 vxl = params(8); vyl = params(9); 

31 vx2 = params(lO); vy2 = params(ll); 

32 vx = (vxl+vx2)/2; vy = (vyl+vy2)/2; 

33 

34 '/, Spacial coordinates in image. 

35 [mz,nz] = size(z); 

36 [x,y] = domain2d([0:mz-l]*delta+vx+l , [0 :nz-l]*delta+vy+l) ; 

37 

38 '/, Stencils 

39 hx = [-1 l]/2/delta; */. X-derivative 

40 hy = [l;0;-l]/2/delta; */. Y-derivative 

41 hxx = [1 -2 1; 2 -4 2; 1 -2 l]/4/delta/delta; 

42 y.hxy = [-1 1;0 0;1 -1] /4/delta/delta; 

43 hxy = [-1 1;1 -1] /delta/delta; 

44 hyy =[12 l;-2 -4 -2 ; 1 2 1] /4/delta/delta; 

45 del2 =[14 1;4 -20 4;1 4 1] /6/delta/delta; 

46 

47 '/, Compute p,q using numerical derivatives of z. 

48 z = hbasis(z,n); 

49 zx = cf ilter2d(z ,hx) ; 

50 zy = cf ilter2d(z ,hy) ; 

51 p = f *zx . /(x. *zx+z) ; 

52 q = f*zy ./(y.*zy+z) ; 

53 

54 '/, Compute reflectance map values. 

55 Rl = rmap(p,q,psl ,qsl) ; 

56 R2 = rmap(p,q,ps2 ,qs2) ; 

57 Rpl = rmapp(p ,q,psl ,qsl) ; 

58 Rp2 = rmapp(p ,q,ps2 ,qs2) ; 

59 Rql = rmapq(p ,q,psl ,qsl) ; 

60 Rq2 = rmapq(p ,q,ps2 ,qs2) ; 

61 

62 '/, Compute numerical derivatives of z. 

63 zxx = cf ilter2d(z ,hxx) ; 

64 zxy = cf ilter2d(z ,hxy ) ; 

65 zyy = cf ilter2d(z ,hyy ) ; 

66 

67 '/, Compute disparity. 

68 d = (f*b/2)./z; 

69 xl = (x+d-vxl-l)/delta+l ; 

70 x2 = (x-d-vx2-l)/delta+l ; 

71 

72 '/, Compute 1st difference in the x direction of image arrays. 

73 '/, Based on the stencil: -1 1 
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74 [ml,nl] = size(El); [m2,n2] = size(E2); 

75 stencil = [-1 l]/delta; 

76 Exl = filter2d(El .stencil, 'resize') ; 

77 Ex2 = filter2d(E2, stencil, 'resize') ; 
78 

79 '/, Determine stereo mapped image derivatives. 

80 Fxl = interpx(Exl ,f loor(xl) ) ; zz = zeros(ml ,nl) ; 

81 out = isnan(Fxl); if any(out(:)), Fxl(out) = zz(out); end 

82 Fx2 = interpx(Ex2 ,f loor(x2) ) ; zz = zeros(m2 ,n2) ; 

83 out = isnan(Fx2) ; if any(out(:)), Fx2(out) = zz(out); end 

84 

85 '/, Determine stereo mapped images F1,F2. Set error to zero where F is NaN 

86 Fl = interpx(El ,xl) ; 

87 out = isnan(Fl); if any(out(:)), Fl(out) = Rl(out); end 

88 F2 = interpx(E2,x2) ; 

89 out = isnan(F2); if any(out(:)), F2(out) = R2(out); end 



ERql = (Fl-Rl).*Rql; 
ERq2 = (F2-R2).*Rq2; 



91 '/, Compute Gradient 

92 G = zeros(mz ,nz) ; 

93 ERpl = (Fl-Rl) .*Rpl; 

94 ERp2 = (F2-R2) .*Rp2; 

95 denx = (x.*zx+z) . ~2 ; 

96 deny = (y . *zy+z) . ~2 ; 

97 G = -(f*b/2)*((Fl-Rl).*Fxl - (F2-R2).*Fx2)./(z."2) + ... 

98 f*((ERpl+ERp2) .*zx./denx + (ERql+ERq2) .*zy . /deny) ; 

99 G(:) = G - cf ilter2d((ERpl+ERp2) . *z . /denx ,-f *hx , 'grad') ; 

100 G(:) = G - cf ilter2d((ERql+ERq2) . *z . /deny ,-f*hy, 'grad') ; 

101 G(:) = G + cf ilter2d(zxx,lambda*hxx, 'grad' ) + ... 

102 cf ilter2d(zxy ,2*lambda*hxy , 'grad' ) + ... 

103 cf ilter2d(zyy ,lambda*hyy , 'grad' ) ; 

104 G = hbasis(G*(l/area) ,n, 'trans ') ; 

B.1.5 hbdfss_cost4.m 



function [J ,FR1 , FR2]=hbdf ss_cost4(v , levels ,params ,E1 ,E2 , lambda, mu) 
HBDFSS_C0ST4 Cost function for depth from shading and stereo 
problem using hierarchical basis functions. Uses true 
perspective projection and dual z maps. 

J=HBDFSS_C0ST4 (V , LEVELS , PARAMS , El , E2 , LAMBDA ,MU) where 
V = [zl ( : ) ;z2( : )] are the optimization variables, LEVELS is the 
number of h-basis levels. El and E2 are the input images, LAMBDA 
is the scalar weighting factor on departure from smoothness and MU 
is the scalar weighting factor on stereo matching. The image 
parameters are PARAMS = [f ,b,z0 ,psl ,qsl ,ps2 ,ps2] . 

J = [Jtot , Jimage , Jsmooth, Jstereo] 

Clay M. Thompson 7-18-91 
Revised to support multi-grid scheme 



2 

3 

4 

5 

6 

7 

8 

9 
10 
11 
12 
13 
14 
15 
16 
17 

18 '/, Camera constants 

19 f = params (1) ; 

20 b = params (2) ; 

21 zO = params(3); 

22 delta = params(12); 

23 area = prod(size(El) )*delta*delta; 

24 

25 '/, Light Source positions 

26 psl = params(4); qsl = params(5); 

27 ps2 = params(6); qs2 = params(7); 

28 

29 '/, Camera coordinate calibration 

30 vxl = params(8); vyl = params(9); 

31 vx2 = params(lO); vy2 = params(ll); 

32 

33 '/, Extract zl & z2 and transform into nodal basis. 

34 [mzl,nzl] = size(El); mzl = mzl+2; nzl = nzl+2; 

35 [mz2,nz2] = size(E2); mz2 = mz2+2 ; nz2 = nz2+2 ; 

36 zl = zeros(mzl ,nzl) ; zl(:) = v(l :mzl*nzl) ; 

37 z2 = zeros (mz2 ,nz2) ; z2(:) = v(mzl*nzl+[l :mz2*nz2] ) ; 

38 zl(:) = hbasis(zl .levels) ; 

39 z2(:) = hbasis(z2 , levels) ; 

40 

41 '/, Spacial coordinates in image. 

42 [xl,yl] = domain2d([0:mzl-l]*delta+vxl + l, [0:nzl-l]*delta+vyl + D ; 

43 [x2,y2] = domain2d([0:mzl-l]*delta+vx2+l, [0:nzl-l]*delta+vy2+l) ; 

44 

45 '/, Stencils 
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46 hx = [-1 l]/2/delta; 

47 hy = [l;0;-l]/2/delta; 

48 del2 = [1 4 1;4 -20 4; 

49 

'/, Compute p,q using numerical derivatives of z. 



'/, X-derivative 
'/. Y-derivative 
1 4 l]/6/delta/delta; 



zxl = f ilter2d(zl ,hx, 
filter2d(z2,hx, 



50 
51 
52 ZX2 

53 
54 
55 
56 
57 
58 
59 
60 
61 



resize 
'resize 



zyl = filter2d(zl ,hy , 'resize ' ) 
zy2 = filter2d(z2 ,hy , 'resize ' ) 



rows = 2:mzl-l; cols = 2:nzl-l 

pi = f*zxl (rows ,:). /(xl (rows , cols) . *zxl (rows ,: )+zl (rows , cols) ) 

ql = f*zyl (: ,cols) . /(yl (rows , cols) . *zyl (: ,cols)+zl (rows ,cols) ) 

rows = 2:mz2-l; cols = 2:nz2-l; 

p2 = f*zx2 (rows ,:). /(x2(rows , cols) . *zx2(rows ,: )+z2(rows , cols) ) 

q2 = f *zy2 ( : ,cols) . /(y2(rows ,cols) . *zy2( : ,cols)+z2(rows ,cols) ) 

'/, Compute reflectance map values. 
Rl = rmap(pl ,ql ,psl ,qsl) ; 



rmap(p2 , q2 ,ps2 , qs2) ; 



62 R2 

63 

64 '/, Compute interpolated images. Set error to zero where zbar is 

65 zbarl = icubic(xl(l,:)',zl' ,(x2+(f*b)./z2) ')'; 

66 out = isnan(zbarl) ; if any(out(:)), zbarl(out) = z2(out); end 

67 zbar2 = icubic(x2(l , : ) ' ,z2 ' , (xl-(f *b) . /zl) ' ) ' 



NaN. 



out = isnan(zbar2) ; if any(out(:)), zbar2(out) 



zl (out) ; end 
'/. sfs 



terml = (0 . 5/area)*sum(sum((El-Rl) . "1 + (E2-R2).~2)) 

term2 = (0 . 5/area)*lambda*sum(sum( ... 

f ilter2d(zl,del2, 'resize') . ~2 + filter2d(z2 ,del2 , 'resize ' ) 

term3 = (0 . 5/area)*mu*sum(sum((zl-zbar2) . ~2 + (zbarl-z2) . ~2) ) ; 
74 '/.disp(sprintf( 'Terms: '/.12.5f '/.12.5f '/.12 .5f ', terml ,term2 ,term3) ) ; 
75 
76 
77 
78 
79 
80 



~2 )) ; '/. Smooth 
'/, Stereo 



J = [terml+term2+term3 , terml ,term2 ,term3] ; 

if nargout>l, 

FR1 = [El;Rl;abs(El-Rl)] ; 

FR2 = [E2;R2;abs(E2-R2)] ; 
end 



B.1.6 hbdfss_grad4.m 



2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
38 
39 
40 



function G=hbdf ss_grad4(v , levels ,params ,E1 ,E2 ,lambda,mu) 
*/,HBDFSS_GRAD4 Gradient function for depth from shading and stereo 

problem using hierarchical basis representation. Uses true 

perspective projection and dual z maps. 

G=HBDFSS_GRAD4 (V , LEVELS , PARAMS , El , E2 , LAMBDA ,MU) where 
V = [zl ( : ) ;z2( : )] are the optimization variables, LEVELS is the 
number of h-basis levels. El and E2 are the input images, LAMBDA 
is the scalar weighting factor on departure from smoothness and MU 
is the scalar weighting factor on stereo matching. The image 
parameters are PARAMS = [f ,b,z0 ,psl ,qsl ,ps2 ,ps2] . 

Clay M. Thompson 7-18-91 

'/, Camera constants 

f = params (1) ; 

b = params (2) ; 

zO = params (3) ; 

delta = params(12); 

area = prod(size(El) )*delta*delta; 

'/. Light Source positions 

psl = params(4); qsl = params(5); 

ps2 = params(6); qs2 = params(7); 

'/, Camera coordinate calibration 
vxl = params(8); vyl = params(9); 
vx2 = params(lO); vy2 = params(ll); 

'/, Extract zl & z2 and transform into nodal basis. 

[mzl,nzl] = size(El); mzl = mzl+2; nzl = nzl+2; 

[mz2,nz2] = size(E2); mz2 = mz2+2 ; nz2 = nz2+2 ; 

zl = zeros(mzl ,nzl) ; zl(:) = v(l :mzl*nzl) ; 

z2 = zeros (mz2 , nz2) ; z2(:) = v(mzl*nzl+[l :mz2*nz2] ) ; 

zl(:) = hbasis(zl , levels) ; 

z2(:) = hbasis(z2 , levels) ; 

'/, Spacial coordinates in image. 

[xl ,yl] = domain2d([0:mzl-l]*delta+vxl + l, [0:nzl-l]*delta+vyl + D ; 

[x2,y2] = domain2d([0:mzl-l]*delta+vx2+l, [0:nzl-l]*delta+vy2+D ; 



42 '/, Stencils 
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'/, X-derivative 
'/. Y-derivative 
1 4 l]/6/delta/delta; '/. Laplacian 



43 hx = [-1 l]/2/delta; 

44 hy = [l;0;-l]/2/delta; 

45 del2 = [1 4 1;4 -20 4; 

46 

47 '/, Compute p,q using numerical derivatives of z. 

48 zxl = filter2d(zl ,hx, 'resize ') ; zyl = filter2d(zl ,hy , 'resize ') ; 

49 zx2 = filter2d(z2 ,hx, 'resize ') ; zy2 = filter2d(z2 ,hy , 'resize ') ; 

50 rows = 2:mzl-l; cols = 2:nzl-l; 

51 pi = f*zxl (rows ,:). /(xl (rows , cols) . *zxl (rows ,: )+zl (rows , cols) ) ; 

52 ql = f*zyl (: ,cols) . /(yl (rows , cols) . *zyl (: ,cols)+zl (rows ,cols) ) ; 

53 rows = 2:mz2-l; cols = 2:nz2-l; 

54 p2 = f*zx2 (rows ,:). /(x2(rows , cols) . *zx2(rows ,: )+z2(rows , cols) ) ; 

55 q2 = f *zy2 ( : ,cols) . /(y2(rows ,cols) . *zy2( : ,cols)+z2(rows ,cols) ) ; 

56 

57 '/, Compute reflectance map values. 

58 Rl = rmap(pl ,ql ,psl ,qsl) ; 

59 Rpl = rmapp(pl ,ql ,psl ,qsl) ; Rql = rmapq(pl ,ql ,psl ,qsl) ; 

60 R2 = rmap(p2 ,q2 ,ps2,qs2) ; 

61 Rp2 = rmapp(p2 ,q2 ,ps2 ,qs2) ; Rq2 = rmapq(p2 ,q2 ,ps2 ,qs2) ; 

62 

63 '/, Compute interpolated images and related terms. Set gradient to zero where zbar is NaN . 

64 d = (f*b)./zl; 

65 zbar2 = icubic(x2(l , : ) ' ,z2 ' , (xl-d) ' ) ' ; 

66 out = f ind(isnan(zbar2) ) ; 

67 zbarx2 = dcubicx(x2(l ,:)' ,z2 ', (xl-d) ') '/delta; 

68 if length (out )>0, zbar2(out) = zl(out); zbarx2(out) = 0*zl(out); end 

69 zbarz2 = dcubicz(x2(l , : ) ' ,z2 ' , (xl-d) ' , (zl-zbar2) ' ) ' ; 
70 

71 d = (f*b)./z2; 

72 zbarl = icubic(xl(l , : ) ' ,zl' , (x2+d) ' ) ' ; 

73 out = f ind(isnan(zbarl) ) ; 

74 zbarxl = dcubicx(xl (1 ,:)' ,zl ', (x2+d) ') '/delta; 

75 if length (out )>0, zbarl(out) = z2(out); zbarxl(out) = 0*z2(out); end 

76 zbarzl = dcubicz(xl (1 , : ) ' ,zl ' , (x2+d) ' , (zbarl-z2) ' ) ' ; 
77 

78 '/, Compute error terms 
79 

80 '/, Now form gradient 

81 rows = 2:mzl-l; cols = 2:nzl-l; 

82 ERp = (El-Rl).*Rpl; ERq = (El-Rl) .*Rql ; 

83 denx = (xl (rows , cols) . *zxl (rows ,: )+zl (rows , cols) ). ~2 ; 

84 deny = (yl (rows , cols) . *zyl (: ,cols)+zl (rows ,cols) ). ~2 ; 

85 dJdzl = zeros (mzl ,nzl ) ; 

86 dJdzKrows , : ) = -filter2d(ERp. *zl (rows , cols) . /denx ,-f *hx) ; 

87 dJdzl( : ,cols) = dJdzl ( : ,cols) - ... 

88 f ilter2d (ERq. *zl (rows ,cols) . /deny ,-f *hy) ; 

89 dJdzl (rows , cols) = dJdzl (rows , cols) + ... 

90 f* (ERp. *zxl (rows, :) ./denx + ERq. *zyl (: ,cols) . /deny) ; 

91 dJdzl(:) = dJdzl + . . . 

92 lambda*filter2d(f ilter2d(zl ,del2, 'resize') ,del2) + ... 

93 mu*(zl-zbar2) .*(l-(f*b)*zbarx2./(zl .~2)) + ... 

94 mu*zbarzl ; 

95 

96 rows = 2:mz2-l; cols = 2:nz2-l; 

97 ERp = (E2-R2) .*Rp2; ERq = (E2-R2) .*Rq2 ; 

98 denx = (x2 (rows , cols) . *zx2(rows ,: )+z2(rows , cols) ). "2 ; 

99 deny = (y2 (rows , cols) . *zy2( : ,cols)+z2(rows , cols) ). "2 ; 

100 dJdz2 = zeros (mz2 ,nz2) ; 

101 dJdz2(rows , : ) = -filter2d(ERp. *z2(rows , cols) . /denx ,-f *hx) ; 

102 dJdz2( : ,cols) = dJdz2 ( : ,cols) - ... 

103 f ilter2d (ERq. *z2 (rows, cols) . /deny ,-f *hy) ; 

104 dJdz2(rows ,cols) = dJdz2 (rows , cols) + ... 

105 f* (ERp. *zx2 (rows, :) ./denx + ERq. *zy2( :, cols) . /deny) ; 

106 dJdz2(:) = dJdz2 + ... 

107 lambda*filter2d(f ilter2d(z2,del2, 'resize') ,del2) - ... 

108 mu*zbarz2 - ... 

109 mu*(zbarl-z2) .*( (f *b)*zbarxl . /(z2 . ~2)+l) ; 
no 

111 '/, Return gradient in h-basis coordinates 

112 dJdzl = hbasis(dJdzl*(l/area) .levels , 'trans' ) ; 

113 dJdz2 = hbasis(dJdz2*(l/area) .levels , 'trans' ) ; 

114 G = [dJdzl(:);dJdz2(:)] ; 

115 

B.1.7 hbdfss cost7.m 

1 function [J ,FR1 , FR2]=hbdf ss_cost (u,n ,params , El ,E2 , lambda) 

2 '/,HBDFSS_C0ST7 Cost function for depth from shading and stereo 

3 '/, problem with heirarchical basis functions. Uses true 

4 '/, perspective projection. 

5 '/. 
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J=HBDFSS_C0ST7(U,N, PARAMS, El, E2, LAMBDA) where U is the disparity, 
and N is the number of H-basis levels . . El and E2 are the input 
images, LAMBDA is the scalar weighting factor on departure from 
smoothness. The image parameters PARAMS = [f ,b ,zO ,psl ,qsl ,ps2 ,ps2] . 

J= [Jtot , Jimage , Jsmooth] . 

Clay M. Thompson 5-23-92 

Uses u = fb/2z - fb/2z_0 to enhance the numerical stability of the algorithm. 

Revised to support multi-grid scheme 



6 

7 

8 

9 
10 
11 
12 
13 
14 
15 
16 

17 '/, Camera constants 

18 f = params (1) ; 

19 b = params (2) ; 

20 zO = params(3); 

21 uO = f*b/2/z0; 

22 delta = params(12); 

23 area = prod(size(u) )*delta*delta; 

24 

25 '/, Light Source positions 

26 psl = params(4); qsl = params(5); 

27 ps2 = params(6); qs2 = params(7); 

28 

29 '/, Camera coordinate calibration 

30 vxl = params(8); vyl = params(9); 

31 vx2 = params(lO); vy2 = params(ll); 

32 vx = (vxl+vx2)/2; vy = (vyl+vy2)/2; 

33 

34 '/, Spacial coordinates in image. 

35 [mz,nz] = size(u); 

36 [x,y] = domain2d([0:mz-l]*delta+vx+l , [0 :nz-l]*delta+vy+l) ; 

37 

38 '/, Stencils 

39 hx = [-1 l]/(2*delta) ; '/. X-derivative 

40 hy = [l;0;-l]/(2*delta) ; '/. Y-derivative 

41 hxx = [1 -2 1; 2 -4 2; 1 -2 l]/(4*delta*delta) ; 

42 '/.hxy = [-1 1;0 0;1 -1] /(4*delta*delta) ; 

43 hxy = [-1 1;1 -l]/(delta*delta) ; 

44 hyy =[12 l;-2 -4 -2 ; 1 2 1] /(4*delta*delta) ; 

45 del2 =[14 1;4 -20 4;1 4 1] /(6*delta*delta) ; 

46 

47 '/, Compute p,q using numerical derivatives of u. 

48 u = hbasis(u,n); 

49 ux = cf ilter2d(u,hx) ; 

50 uy = cf ilter2d(u,hy) ; 

51 p = f *ux./(x.*ux-u-u0) ; 

52 q = f *uy ./(y .*uy-u-u0) ; 

53 

54 '/, Compute reflectance map values. 

55 Rl = rmap(p,q,psl ,qsl) ; 

56 R2 = rmap(p,q,ps2 ,qs2) ; 

57 

58 '/, Compute numerical derivatives of u. 

59 uxx = cf ilter2d(u ,hxx) ; 

60 uxy = cf ilter2d(u ,hxy ) ; 

61 uyy = cf ilter2d(u ,hyy ) ; 

62 

63 '/, Determine stereo mapped images F1,F2. Set error to zero where F is NaN 

64 Fl = interpx(El, (x+u+uO-vxl-l)/delta+l) ; 

65 out = isnan(Fl); if any(out ( : )) , Fl(out) = Rl(out); end, 

66 F2 = interpx(E2, (x-u-u0-vx2-l)/delta+l) ; 

67 out = isnan(F2) ; if any(out(:)), F2(out) = R2(out); end, 

68 

69 '/.image([Fl,F2;Rl,R2;abs(Fl-Rl),abs(F2-R2)] , [0 1]) 
70 

71 terml = (0 . 5/area)*sum(sum((Fl-Rl) . ~2 + (F2-R2) . ~2) ) ; 

72 term2 = (0 . 5/area)*lambda*(sum(sum(uxx . ~2) ) + 2*sum(sum(uxy . ~2) ) + sum(sum(uyy . "2) ) ) ; 

73 */,term2 = (0 .5/area)*lambda*sum(sum(delsqz . "2) ) ; 

74 '/.disp(sprintf ( 'Terms: '/.12.5f '/.12 . 5f ', terml ,term2) ) ; 

75 J = [terml+term2 , terml ,term2] ; 
76 

77 if nargout>l, 

78 FR1 = [El;Fl;Rl;abs(Fl-RD] ; 

79 FR2 = [E2;F2;R2;abs(F2-R2)] ; 

80 end 



B.1.8 hbdfss grad7.m 

l function G=hbdf ss_grad(u , n, params ,E1 ,E2 , lambda) 
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12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
38 
39 
40 
41 
42 
43 
44 



*/,HBDFSS_GRAD7 Gradient function for depth from shading and stereo 

problem with heirarchical basis functions. Uses true perspective 
projection . 

G=HBDFSS_GRAD7(U,N, PARAMS, El, E2 , LAMBDA) where U is the disparity 
map, and N is the number of H-basis levels. El and E2 are the input 
images, LAMBDA is the scalar weighting factor on departure from 
smoothness. The image parameters PARAMS = [f ,b ,z0 ,psl ,qsl ,ps2 ,ps2] ; 

Clay M. Thompson 5-23-92 

Uses u = fb/2z - fb/2z_0 to enhance the numerical stability of the algorithm. 

Revised to support multi-grid algorithm 

'/, Camera constants 

f = params (1) ; 

b = params (2) ; 

zO = params (3) ; 

uO = f*b/2/z0; 

delta = params(12); 

area = prod(size(u) )*delta*delta; 

'/, Light Source positions 

psl = params(4); qsl = params(5); 

ps2 = params(6); qs2 = params(7); 

'/, Camera coordinate calibration 
vxl = params(8); vyl = params(9); 
vx2 = params(lO); vy2 = params(ll); 
vx = (vxl+vx2)/2; vy = (vyl+vy2)/2; 

'/, Spacial coordinates in image. 

[mz,nz] = size(u); 

[x,y] = domain2d([0:mz-l]*delta+vx+l , [0 :nz-l]*delta+vy+l) ; 

'/, Stencils 

hx = [-1 l]/2/delta 

hy = [l;0;-l]/2/delta 

hxx = [1-2 1; 2-4 2 

'/.hxy = [-1 1;0 

hxy = [-1 1;1 

hyy = [1 2 1; 



'/, X-derivative 
'/. Y-derivative 
1 -2 l]/4/delta/delta; 
1 -l]/4/delta/delta; 
-1] /delta/delta; 
2 -4 -2;1 2 l]/4/delta/delta; 



del2 =[14 1;4 -20 4;1 4 1] /6/delta/delta; 



45 '/, Compute p,q using numerical derivatives of z. 
u = hbasis (u,n) ; 
uu = u + uO ; 
ux = cf ilter2d(u ,hx) ; 
uy = cf ilter2d(u ,hy) ; 
denx = (x.*ux-uu); 
deny = (y.*uy-uu); 
p = f *ux . /denx ; 
q = f *uy . /deny ; 

'/, Squared denominators for gradient terms below 
denx(:) = denx.+denx; 
deny(:) = deny.+deny; 

'/, Compute reflectance map values. 

Rl = rmap(p,q,psl ,qsl) ; 

R2 = rmap(p,q,ps2 ,qs2) ; 

Rpl = rmapp(p ,q,psl ,qsl) ; 

Rp2 = rmapp(p ,q,ps2 ,qs2) ; 

Rql = rmapq(p ,q,psl ,qsl) ; 

Rq2 = rmapq(p ,q,ps2 ,qs2) ; 

'/, Compute numerical derivatives of z. 
uxx = cf ilter2d(u ,hxx) ; 
uxy = cf ilter2d(u ,hxy ) ; 
uyy = cf ilter2d(u ,hyy ) ; 



'/, Compute image projections, 
xl = (x+uu-vxl-l)/delta+l ; 
x2 = (x-uu-vx2-l)/delta+l ; 

'/, Compute 1st difference in the x direction of image arrays. 

'/, Based on the stencil: -1 1 

[ml,nl] = size(El); [m2,n2] = size(E2); 

stencil = [-1 l]/delta; 

Exl = filter2d(El .stencil, 'resize') ; 

Ex2 = filter2d(E2 , stencil , 'resize ') ; 
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82 

83 '/, Determine stereo mapped image derivatives. 

84 Fxl = interpx(Exl ,f loor(xl) ) ; zz = zeros(ml ,nl) ; 

85 out = isnan(Fxl); if any(out(:)), Fxl (out) = zz(out); end 

86 Fx2 = interpx(Ex2 ,f loor(x2) ) ; zz = zeros(m2 ,n2) ; 

87 out = isnan(Fx2) ; if any(out(:)), Fx2(out) = zz(out); end 

88 

89 '/, Determine stereo mapped images F1,F2. Set error to zero where F is NaN 

90 Fl = interpx(El ,xl) ; 

91 out = isnan(Fl); if any(out(:)), Fl(out) = Rl(out); end 

92 F2 = interpx(E2,x2) ; 

93 out = isnan(F2); if any(out(:)), F2(out) = R2(out); end 

94 

95 '/, Compute Gradient 

96 G = zeros(mz ,nz) ; 

97 ERpl = (Fl-Rl) .*Rpl; ERql = (Fl-Rl) . *Rql ; 

98 ERp2 = (F2-R2) .*Rp2; ERq2 = (F2-R2) . *Rq2 ; 

99 G = ((Fl-Rl). *Fxl - (F2-R2).*Fx2) ... 
- f*((ERpl+ERp2) .*ux./denx + (ERql+ERq2) . *uy . /deny) ; 

G(:) = G + cf ilter2d((ERpl+ERp2) . *uu . /denx ,-f *hx, 'grad') ; 
G(:) = G + cf ilter2d((ERql+ERq2) . *uu . /deny ,-f*hy, 'grad') ; 
G(:) = G + cf ilter2d(uxx,lambda*hxx, 'grad') + ... 

cf ilter2d(uxy ,2*lambda*hxy , 'grad ' ) + ... 

cf ilter2d(uyy ,lambda*hyy , 'grad' ) ; 



100 
101 
102 
103 
104 
105 
106 

107 G = hbasis(G*(l/area) ,n, 'trans ') ; 

B.1.9 hbdfss_cost8b.m 



6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 



function [J ,FR1 , FR2]=hbdf ss_cost (v,n ,params ,E1 , E2 , lambda, mu) 
HBDFSS_C0ST8B Cost function for depth from shading and stereo 
problem with heirarchical basis functions. Uses true 
perspective projection and varying albedo 

J=HBDFSS_C0ST8B(V,N, PARAMS, El, E2 , LAMBDA, MU) where V=[Z,RH0] are 
the optimization variables (height and albedo), N is the number of 
H-basis levels.. El and E2 are the input images, LAMBDA and MU are 
scalar weighting factors on departure from smoothness for height 
and albedo, respectively. The image parameters PARAMS = 
[f ,b,z0,psl ,qsl ,ps2,ps2] . 

J= [Jtot , J image , Jsmooth , Jalbedo] . 
Clay M. Thompson 6-28-92 



Camera constants 

18 f = params (1) ; 

19 b = params (2) ; 

20 zO = params(3); 

21 gamma = f*b/2/z0; 

22 delta = params(12); 

23 area = prod(size(El) )*delta*delta; 

24 

25 '/, Light Source positions 

26 psl = params(4); qsl = params(5); 

27 ps2 = params(6); qs2 = params(7); 

28 

29 '/, Camera coordinate calibration 

30 vxl = params(8); vyl = params(9); 

31 vx2 = params(lO); vy2 = params(ll); 

32 vx = (vxl+vx2)/2; vy = (vyl+vy2)/2; 

33 

34 '/, Spacial coordinates in image. 

35 [mz,nz] = size (El); 

36 z = v(:,l:nz); rho = v( : ,nz+[l :nz] ) ; 

37 [x,y] = domain2d([0:mz-l]*delta+vx+l , [0 :nz-l]*delta+vy+l) ; 

38 

39 '/, Stencils 

40 hx = [-1 l]/(2*delta) ; '/. X-derivative 

41 hy = [l;0;-l]/(2*delta) ; '/. Y-derivative 

42 hxx = [1-2 1; 2-42; 1 -2 l]/(4*delta*delta) ; 

43 '/.hxy = [-1 1;0 0;1 -1] /(4*delta*delta) ; 

44 hxy = [-1 1;1 -l]/(delta*delta) ; 

45 hyy =[12 l;-2 -4 -2 ; 1 2 1] /(4*delta*delta) ; 

46 del2 =[14 1;4 -20 4;1 4 1] /(6*delta*delta) ; 

47 hr = [0 0;0 1 0;0 0] - ones(3,3)/9; 

48 

49 '/, Compute p,q using numerical derivatives of z. 
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50 z = hbasis(z,n); rho = hbasis(rho ,n) ; 

51 zx = cf ilter2d(z ,hx) ; 

52 zy = cf ilter2d(z ,hy) ; 

53 p = f *zx./(x.*zx+z) ; 

54 q = f *zy . /(y . *zy+z) ; 

55 

56 '/, Compute reflectance map values. 

57 Rl = rho . *rmap(p ,q,psl ,qsl) ; 

58 R2 = rho . *rmap(p ,q,ps2 ,qs2) ; 

59 

60 '/, Compute numerical derivatives of z. 

61 zxx = cf ilter2d(z ,hxx) ; 

62 zxy = cf ilter2d(z ,hxy) ; 

63 zyy = cf ilter2d(z ,hyy) ; 

64 

65 '/, Compute disparity. 

66 d = (f*b/2)./z; 

67 

68 '/, Determine stereo mapped images F1,F2. Set error to zero where F is NaN 

69 Fl = interpx(El, (x+d-vxl-l)/delta+l) ; 

70 out = isnan(Fl); if any(out(:)), Fl(out) = Rl(out); end, 

71 F2 = interpx(E2, (x-d-vx2-l)/delta+l) ; 

72 out = isnan(F2); if any(out(:)), F2(out) = R2(out); end, 
73 

74 '/.image ( [Fl, F2 ; Rl, R2] , [0 1],5) 
75 

76 terml = (0 . 5/area)*sum(sum((Fl-Rl) . ~2 + (F2-R2) . ~2) ) ; 

77 term2 = (0 . 5/area)*lambda*(sum(sum(zxx . ~2) ) + 2*sum(sum(zxy . ~2) ) + sum(sum(zyy . "2) ) ) ; 

78 term3 = (0 . 5/area)*mu*sum(sum(cf ilter2d(rho ,hr) . "2) ) ; 

79 */,term2 = (0 .5/area)*lambda*sum(sum(delsqz . "2) ) ; 

80 '/.disp(sprintf ( 'Terms: '/.12.5f '/.12.5f '/.12 .5f ', terml ,term2 ,term3) ) ; 

81 J = [terml+term2+term3 , terml ,term2 ,term3] ; 

82 

83 if nargout>l, 

84 FR1 = [El;Fl;Rl;abs(F2-R2)] ; 

85 FR2 = [E2;F2;R2;abs(F2-R2)] ; 

86 end 



B.1.10 hbdfss_grad8b.m 



6 

7 

8 

9 

10 

11 

12 

13 

14 

15 



function G=hbdf ss_grad(v ,n,params ,E1 ,E2 ,lambda,mu) 

*/,HBDFSS_GRAD8B Gradient function for depth from shading and stereo 

problem with heirarchical basis functions. Uses true perspective 

projection and varying albedo. 

G=HBDFSS_GRAD8B(V,N,PARAMS, El, E2, LAMBDA, MU) where V=[Z,RH0] are 
the optimization variables (height and albedo), N is the number of 
H-basis levels.. El and E2 are the input images, LAMBDA and MU are 
scalar weighting factors on departure from smoothness for height 
and albedo, respectively. The image parameters PARAMS = 
[f ,b,z0,psl ,qsl ,ps2,ps2] . 



Clay M. Thompson 6-28-92 



Camera constants 

16 f = params(l) ; 

17 b = params (2) ; 

18 zO = params(3); 

19 gamma = f*b/2/z0; 

20 delta = params(12); 

21 area = prod(size(El) )*delta*delta; 

22 

23 '/, Light Source positions 

24 psl = params(4); qsl = params(5); 

25 ps2 = params(6); qs2 = params(7); 

26 

27 '/, Camera coordinate calibration 

28 vxl = params(8); vyl = params(9); 

29 vx2 = params(lO); vy2 = params(ll); 

30 vx = (vxl+vx2)/2; vy = (vyl+vy2)/2; 

31 

32 '/, Spacial coordinates in image. 

33 [mz,nz] = size (El); 

34 z = v(:,l:nz); rho = v( : ,nz+[l :nz] ) ; 

35 [x,y] = domain2d([0:mz-l]*delta+vx+l , [0 :nz-l]*delta+vy+l) ; 

36 

37 '/, Stencils 

38 hx = [-1 l]/2/delta; '/. X-derivative 

39 hy = [l;0;-l]/2/delta; '/. Y-derivative 
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40 hxx = [1-2 1; 2-42; 1 -2 l]/4/delta/delta; 

41 '/.hxy = [-1 1;0 0;1 -1] /4/delta/delta; 

42 hxy = [-1 1;1 -1] /delta/delta; 

43 hyy = [1 2 l;-2 -4 -2 ; 1 2 1] /4/delta/delta; 

44 del2 =[14 1;4 -20 4;1 4 1] /6/delta/delta; 

45 hr = [0 0;0 1 0;0 0] - ones(3,3)/9; 

46 

47 '/, Compute p,q using numerical derivatives of z. 

48 z = hbasis(z,n); rho = hbasis(rho ,n) ; 

49 zx = cf ilter2d(z ,hx) ; 

50 zy = cf ilter2d(z ,hy) ; 

51 p = f *zx . /(x. *zx+z) ; 

52 q = f *zy ./(y .*zy+z) ; 

53 

54 '/, Compute reflectance map values. 

55 RR1 = rmap(p,q,psl ,qsl) ; Rl = rho . *RR1 ; 

56 RR2 = rmap(p,q,ps2,qs2) ; R2 = rho . *RR2 ; 

57 Rpl = rho . *rmapp(p,q,psl ,qsl) ; 

58 Rp2 = rho . *rmapp(p,q,ps2 ,qs2) ; 

59 Rql = rho . *rmapq(p,q,psl ,qsl) ; 

60 Rq2 = rho . *rmapq(p,q,ps2 ,qs2) ; 

61 

62 '/, Compute numerical derivatives of z. 

63 zxx = cf ilter2d(z ,hxx) ; 

64 zxy = cf ilter2d(z ,hxy ) ; 

65 zyy = cf ilter2d(z ,hyy ) ; 

66 

67 '/, Compute disparity. 

68 d = (f*b/2)./z; 

69 xl = (x+d-vxl-l)/delta+l ; 

70 x2 = (x-d-vx2-l)/delta+l ; 
71 

72 '/, Compute 1st difference in the x direction of image arrays. 

73 '/, Based on the stencil: -1 1 

74 [ml,nl] = size(El); [m2,n2] = size(E2); 

75 stencil = [-1 l]/delta; 

76 Exl = filter2d(El .stencil, 'resize') ; 

77 Ex2 = filter2d(E2, stencil, 'resize') ; 
78 

79 '/, Determine stereo mapped image derivatives. 

80 Fxl = interpx(Exl ,f loor(xl) ) ; zz = zeros(ml ,nl) ; 

81 out = isnan(Fxl); if any(out(:)), Fxl(out) = zz(out); end 

82 Fx2 = interpx(Ex2 ,f loor(x2) ) ; zz = zeros(m2 ,n2) ; 

83 out = isnan(Fx2) ; if any(out(:)), Fx2(out) = zz(out); end 

84 

85 '/, Determine stereo mapped images F1,F2. Set error to zero where F is NaN 

86 Fl = interpx(El ,xl) ; 

87 out = isnan(Fl); if any(out(:)), Fl(out) = Rl(out); end 

88 F2 = interpx(E2,x2) ; 

89 out = isnan(F2) ; if any(out(:)), F2(out) = R2(out); end 

90 

91 '/, Compute Gradient 

92 FR1 = (Fl-Rl); FR2 = (F2-R2) ; 

93 ERpl = FRl.*Rpl; ERql = FR1 . *Rql ; 

94 ERp2 = FR2 . *Rp2 ; ERq2 = FR2 . *Rq2 ; 

95 denx = (x.*zx+z) . ~2 ; 

96 deny = (y . *zy+z) . ~2 ; 

97 Gz = -(f*b/2)*(FRl.*Fxl - FR2 . *Fx2) . /(z . ~2) + ... 

98 f*((ERpl+ERp2) .*zx./denx + (ERql+ERq2) .*zy . /deny) ; 

99 Gz(:) = Gz - cfilter2d((ERpl+ERp2) .*z./denx,-f*hx, 'grad') ; 

100 Gz(:) = Gz - cfilter2d((ERql+ERq2) . *z . /deny, -f *hy , 'grad') ; 

101 Gz(:) = Gz + cf ilter2d (zxx, lambda*hxx, 'grad' ) + ... 

102 cf ilter2d(zxy ,2*lambda*hxy , 'grad' ) + ... 

103 cf ilter2d(zyy ,lambda*hyy , 'grad' ) ; 

104 

105 Grho = cfilter2d(cf ilter2d(rho ,hr) ,mu*hr, 'grad') - FR1 . *RR1 - FR2 . *RR2 ; 

106 

107 Gz = hbasis(Gz*(l/area) ,n, 'trans' ) ; 

108 Grho = hbasis(Grho*(l/area) ,n, 'trans ') ; 

109 G = [Gz.Grho] ; 

no 

B.l.ll hbsfs_cost.m 

1 function [J ,FR]=hbsf s_cost (z , levels ,params ,E , lambda) 

2 '/.HBSFS Cost function for shape from shading using hierarchical 

3 '/, basis functions. Based on orthographic projection. 

4 '/. 

5 '/. J=HBSFS_C0ST(Z, LEVELS, PARAMS.E, LAMBDA) where Z is the depth map, 
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6 '/, LEVELS is the number of h-basis levels. E is the input image, and 

7 '/, LAMBDA is the scalar weighting factor on departure from 

8 '/, smoothness. The image parameters are PARAMS = [f ,b,zO,ps ,qs] . 

9 

10 */. Clay M. Thompson 12-16-91 
11 

12 '/, Camera constants 

13 f = params (1) ; 

14 b = params (2) ; 

15 zO = params(3); 

16 delta = params(12); 

17 area = prod(size(z) )*delta*delta; 

18 

19 '/, Light Source positions 

20 ps = params(4); qs = params(5); 

21 

22 '/, Extract zl & z2 and transform into nodal basis. 

23 [mz,nz] = size(E); mz = mz+2 ; nz = nz+2 ; 

24 z = hbasis (z , levels) ; 

25 

26 '/, Stencils 

27 hx = (f/z0)*[-l l]/2/delta; */. X-derivative 

28 hy = (f/z0)*[l ;0;-l]/2/delta; */. Y-derivative 

29 del2 = [1 4 1;4 -20 4;1 4 1] /6/delta/delta; */. Laplacian 

30 

31 '/, Compute p,q using numerical derivatives of z. 

32 rows = 2:mz-l; cols = 2:nz-l; 

33 p = f ilter2d(z ,hx , 'resize ') ; p = p(rows,:); 

34 q = f ilter2d(z ,hy , 'resize ') ; q = q(:,cols); 

35 

36 '/, Compute reflectance map values. 

37 R = rmap(p ,q,ps ,qs) ; 

38 

39 terml = (0 . 5/area)*sum(sum((E-R) . "2) ) ; '/. sfs 

40 term2 = (0 . 5/area)*lambda*sum(sum( f ilter2d(z ,del2 , 'resize '). "2 )); '/, Smoothness 

41 */.disp(sprintf( 'Terms: */.12 .5f */.12 . 5f ' , terml ,term2) ) ; 

42 J = [terml+term2 ,terml ,term2] ; 

43 

44 if nargout>l, 

45 FR = [E;R;abs(E-R>] ; 

46 end 



B.1.12 hbsfs_grad.m 



function G=hbsf s_grad(z , levels , params ,E , lambda) 
'/,HBDFSS_GRAD4 Gradient function for shape from shading using 

hierarchical basis representation. Based on orthographic 

projection . 

G=HBSFS_GRAD(Z, LEVELS, PARAMS, E, LAMBDA) where Z is the depth map, 
LEVELS is the number of h-basis levels. E is the input image, and 
LAMBDA is the scalar weighting factor on departure from 
smoothness. The image parameters are PARAMS = [f ,b,z0,ps ,qs] . 



Clay M. Thompson 12-16-91 

Revised to support multi-grid scheme. 



10 

n 

12 
13 

14 '/, Camera constants 

15 f = params (1) ; 

16 b = params (2) ; 

17 zO = params(3); 

18 delta = params(12); 

19 area = prod(size(z) )*delta*delta; 

20 

21 '/, Light Source positions 

22 ps = params(4); qs = params(5); 

23 

24 '/, Extract zl & z2 and transform into nodal basis. 

25 [mz,nz] = size(E); mz = mz+2; nz = nz+2; 

26 z = hbasis (z , levels) ; 

27 

28 '/, Stencils 

29 hx = (f/z0)*[-l l]/2/delta; V. X-derivative 

30 hy = (f/z0)*[l ;0;-l]/2/delta; V. Y-derivative 

31 del2 = [1 4 1;4 -20 4;1 4 1] /6/delta/delta; */. Laplacian 

32 

33 '/, Compute p,q using numerical derivatives of z. 

34 rows = 2:mz-l; cols = 2:nz-l; 

35 p = f ilter2d(z ,hx , 'resize ') ; p = p(rows,:); 

36 q = f ilter2d(z ,hy , 'resize ') ; q = q(:,cols); 
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37 del2z = filter2d(z,del2, 'resize') ; 

38 

39 '/, Compute reflectance map values. 

40 R = rmap(p ,q,ps ,qs) ; 

41 Rp = rmapp(p,q,ps ,qs) ; Rq = rmapq(p ,q,ps ,qs) ; 

42 

43 '/, Compute error terms 

44 

45 '/, Now form gradient 

46 G = zeros(mz ,nz) ; 

47 GCrows,:) = f ilter2d( (E-R) . *Rp ,hx) ; 

48 G(:, cols) = G(:, cols) + f ilter2d((E-R) .*Rq,hy) ; 

49 GO) = G + lambda*f ilter2d(del2z ,del2) ; 

50 

51 '/, Return gradient in h-basis coordinates 

52 G = hbasis(G*(l/area) .levels , 'trans ') ; 



B.2 Support routines. 

B.2.1 rmap.m 

l function R=rmap(p ,q,ps ,qs) 



'/.RMAP Reflectance map calculation 

'/, R=RMAP(P,Q) computes the reflectance map image of the surface with 

'/. the gradients P and Q (in the x and y direction respectively). 

'/, P and Q are matrices that contain the gradients over a rectangular 

*/. grid. 

% 

*/. R=RMAP(P,Q,Ps,Qs) uses the light source direction (Ps,Qs). 

% 

'/, Currently implements: Lambertian reflectance 



12 if nargin==2 , 

13 ps = . 1 ; qs = . 1 ; '/, Light source direction 

14 end 

15 

16 [n ,m] = size(p) ; 

17 R = max( (l+ps*p+qs*q) ./ sqrt (1+p. *p+q.*q) ./ sqrt (l+ps*ps+qs*qs) 

18 ,zeros(n ,m) ) ; 



B.2. 2 rmapp.m 



1 function Rp=rmapp(p ,q,ps ,qs) 

2 '/.RMAPP Reflectance map partial derivative calculation 

3 '/, Rp=RMAPP(P,Q) computes the reflectance map X partial derivative of 

4 '/, the surface with gradients P and Q (in the x and y direction 

5 '/, respectively) . P and Q are matrices that contain the gradients 

6 '/, over a rectangular grid. 

7 '/. 

8 '/. Rp=RMAPP(P,Q,Ps,Qs) uses the light source direction (Ps,Qs). 

9 '/. 

10 '/, Currently implements: Lambertian reflectance 
11 

12 if nargin==2 , 

13 ps = . 1 ; qs = . 1 ; '/, Light source direction. 

14 end 

15 

16 [n ,m] = size(p) ; 

17 d = ones(n,m) + p.*p + q.*q; 

18 e = ones(n,m) + ps*p + qs*q; 

19 Rp = (ps - e .* p ./ d ) ./ sqrt(d) ./ sqrt (l+ps*ps+qs*qs) ; 

20 

21 ndx = f ind(rmap(p ,q,ps ,qs)==0) ; 

22 if length(ndx)X) , Rp(ndx) = zeros(length(ndx) ,1) ; end 

23 



B.2. 3 rmapq.m 



1 function Rq=rmapq(p ,q,ps ,qs) 

2 '/.RMAPQ Reflectance map partial derivative calculation 

3 '/, Rq=RMAPQ (P,Q) computes the reflectance map Y partial derivative of 

4 '/, the surface with gradients P and Q (in the x and y direction 

5 '/, respectively) . P and Q are matrices that contain the gradients 

6 '/, over a rectangular grid. 
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r '/. 

8 '/. Rq=RMAPQ(P,Q,Ps,Qs) uses the light source direction (Ps,Qs). 

9 '/. 

10 '/, Currently implements: Lambertian reflectance 
11 

12 if nargin==2 , 

13 ps = . 1 ; qs = . 1 ; '/, Light source direction. 

14 end 

15 

16 [n ,m] = size(p) ; 

17 d = ones(n,m) + p.*p + q.*q; 

18 e = ones(n,m) + ps*p + qs*q; 

19 Rq = (qs - e . * q ./ d) ./ sqrt(d) ./ sqrt (l+ps*ps+qs*qs) ; 

20 

21 ndx = f ind(rmap(p ,q,ps ,qs)==0) ; 

22 if length(ndx)>0 , Rq(ndx) = zerosdength(ndx) ,1) ; end 

23 



B.2.4 conjgrad.m 



1 function [x .OPTIONS , J ,history]=conjgrad(FUN,x, OPTIONS , GRAD , . 

2 P1,P2,P3,P4,P5,P6,P7,P8,P9,P10) 
'/.CONJGRAD Conjugate-gradient optimization. 



X = CONJGRAD ('FUN' ,X0, OPTIONS, 'GRAD') finds the minimum of the 
function 'FUN' with gradient 'GRAD' using a conjugate-gradient 
optimization algorithm. XO is an initial condition. OPTIONS is a 
vector that contains optional information for the optimizer (see 
F0PTI0NS). 'FUN' and 'GRAD' are strings that contain the names of 
the cost function and gradient M-f iles , respectively. 'FUN' should 
return a scalar function value, f=FUN(x) . 'GRAD' should return the 
gradient vector (df /dx) , g = GRAD(x) . 



Up to ten parameters can be passed to 'FUN' and 'GRAD' using X = 
CONJGRAD ('fun'.XO, OPTIONS, 'grad' ,P1,P2,. . .) so f =FUN(x ,P1 ,P2 , . . .) 
and g=GRAD(x,Pl,P2,. . .). 

Note: FUNCTION value is returned in 0PTI0NS(8) and gradient value 
is returned in 0PTI0NS(15). 0PTI0NS(19)-History save rate. 

'/. Clay M. Thompson 2-4-91 



12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 

23 tol = 1 . e-6 ; '/, Minimum allowed alpha 

24 

25 error(nargchk(4, 14,nargin) ) ; 

26 

27 '/, Form call strings. 

28 params = [] ; 

29 for n=5:nargin 

30 params = [params ,' ,P ', int2str(n-4)] ; 

31 end 

32 if ~any(FUN<48), fcall = [FUN, ' (x ', params ,')'] ; else fcall = FUN; end 

33 if ~any(GRAD<48) , gcall = [GRAD ,' (x ', params ,')'] ; else gcall = GRAD; end 

34 if ~any(FUN<48), 

35 LPARAM = [] ; 

36 for n=5:nargin, 

37 LPARAM = [LPARAM, ' ,P' ,int2str(n-2)] ; 

38 end 

39 SEARCH = [' [eval("x=Pl+x*P2; "),', FUN, ' (x ' , LPARAM, ')]'] ; 

40 si S 6 

41 SEARCH = [' [eval("x=Pl+x*P2 ;"),', FUN, ']'] ; 

42 end 

43 linecall = ['lsearch (SEARCH ,0,2*alpha,f ,GRAD ,x ,pk ' , params , ' ) '] ; 

44 

45 nvars = length (x( : ) ) ; 

46 [mx,nx] = size(x); 

47 

48 '/, Initialize parameters 

49 beta = 0; 

50 f = eval (fcall); '/, Function value 

51 fold = f ; 

52 gk = eval(gcall); '/, Gradient 

53 pk = zeros (mx ,nx) ; '/, Search direction 

54 gknorm = norm(gk(:)); 

55 OPTIONS = foptions(OPTIONS); 

56 OPTIONS(IO) = OPTIONS(IO) + 1; 

57 OPTIONS(ll) = OPTIONS(ll) + 1; 

58 if OPTIONS (18)==0, alpha=.01; else, alpha = 2*0PTI0NS(18) ; end '/. Initial guess 

59 

60 if size(gk)~=size(x) , error('The size of the gradient and x don''t match.'); end 
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61 if OPTIONS(1)>0, 

62 disp('Fun. Evals Value Stepsize Gradient'); 

63 disp([sprintf("/.5.0f '/.12.3g '/.12.3g ' ,OPTIONS(10) ,f(l) .OPTIONS (18) ) , ... 

64 sprintf C/.12.3g ' ,gknorm)] ) ; 

65 end 

66 if OPTIONS (19) >0 & nargout>3, 

67 history = [OPTIONS(IO) ,f ,OPTIONS(10) ,gknorm] ; 

68 end 

69 

70 while OPTIONS(10)<OPTIONS(14), 

71 '/, Test for convergence . 

72 if gknorm<0PTI0NS(2), 

73 0PTI0NS(8) = f(l); 0PTI0NS(15) = gknorm; 

74 disp( 'Gradient criteria met.'), break 

75 end 

76 if OPTIONS (7)==0, 

77 pk = -gk + beta*pk; '/, Compute search direction. 

78 else 

79 pk = -gk; 

80 end 

81 

82 '/, Do an inexact line search to determine, alpha. 

83 GRAD = gk(:) '*pk(:); 

84 if GRAD>0, 

85 pk = -gk; GRAD = -gk( : ) '*gk( : ) ; 

86 if OPTIONS(1)>0, disp( 'Redirect search.'), end, 

87 end 

88 [alpha, f ,n, how] = eval (linecall) ; '/, Line search 

89 

90 OPTIONS(IO) = OPTIONS(IO) + n; 

91 0PTI0NS(18) = alpha; 

92 

93 x = x + alpha*pk; 

94 

95 if alpha>tol, 

96 '/, gkold = gk ; '/, Needed only with Polak-Ribiere Method 

97 gk = eval(gcall); OPTIONS(ll) = OPTIONS(ll) + 1; 

98 

99 '/.beta = ( (gk-gkold) '*gk) /gknorm; '/, Polak-Ribiere Method 

100 '/.gknorm = norm(gk(:)); 

101 

102 beta = norm(gk( : ) )/gknorm; '/, Fletcher-Reeves Method 

103 gknorm = beta+gknorm; 

104 else 

105 beta = 0; 

106 alpha = tol; 

107 gknorm = norm(gk(:)); 

108 if 0PTI0NS(l)>O, dispCReset alpha.'), end 

109 end 
no 

111 if 0PTI0NS(l)>O, 

112 disp( [sprintf ("/.5. Of '/.12.3g '/.12.3g ' ,0PTI0NS(lO) ,f (1) ,0PTI0NS(18) ) , ... 

113 sprintf ( ''/,12 .3g ', gknorm) , how] ) ; 

114 end 

115 

116 if OPTIONS(19)>0 & nargout>3, 

117 if 0PTI0NS(10)>=history(length(history(: ,1 ) ) , 1)+0PTI0NS(19) , 

118 history = [history ; [OPTIONS(IO) ,f ,0PTI0NS(18) , gknorm]] ; 

119 end 

120 end 

121 

122 '/, Test for convergence on relative change in F. 

123 if abs((fold(l)-f (1 ) )/(f old(l)+eps) )<0PTI0NS (3) , 

124 0PTI0NS(8) = f(l); 0PTI0NS(15) = gknorm; 

125 disp( 'Relative function change criteria met.'), break 

126 end 

127 fold = f ; 

128 end 

129 0PTI0NS(8) = f (1) ; 

130 0PTI0NS(15) = gknorm; 

131 J = f ; 

B.2.5 lsearch.m 

1 function [alpha, f ,n ,how] = lsearch(FUN ,x0 ,alphal ,f 0,g0,Pl ,P2 ,P3 ,P4 ,P5 ,P6 ,P7 ,P8 ,P9 ,P10 ,P11 ,P12) 

2 '/.LSEARCH Inexact line search. 

3 V. 

4 */. [ALPHA, F,N] = LSEARCH('FUN' ,X , ALPHA ,F ,G) performs an inexact 

5 '/, line search of the 1-D function 'FUN'. 'FUN' is a string variable 

6 '/, that defines the name of an M-file function to be minimized. The 

7 '/, function should return a scalar value, f=FUN(x) . X is the initial 

8 '/, starting point. F and G are the value and gradient of the function 
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10 

11 

12 
13 
14 



'/, at the initial point X. G must be negative (i.e., downhill). ALPHA 

'/, is an initial guess for the step length (i.e. XI = X + ALPHA). 

'/, LSEARCH returns the step length ALPHA, the function value F, and 

'/, the number of function evaluations N. 

'/. 

'/, Up to ten additianal parameters can be passed to 'fun' using 

'/. [ALPHA, F,N] = LSEARCH (> FUN' ,X, ALPHA, F,G, PI, P2,P3, .. .) in which case 

'/. the function is called using FUN(x ,P1 ,P2 ,P3 , . . . ) . 



16 
17 

18 '/. Clay M. Thompson 4-11-91 

19 

20 showSteps = 0; 

21 showPlot = 0; 

22 

23 tol = le-7; '/, Minimum tolerable value for alpha 

24 

25 '/, Check for valid initial condition. 

26 if g0>0, error( ' Initial gradient is not downhill.'); end 

27 

28 '/, Form call string 

29 params = ' (x' ; 

30 for n=6 margin 

31 params = [params ,' ,P ', int2str(n-5)] ; 

32 end 

33 if ~any(FUN<48), fcall = [FUN, params , ' ) '] ; else fcall = FUN; end 

34 

35 how = ' ' ; 

36 n = 0; 

37 

38 if showSteps, disp(' Alpha Function Value'), end 

39 if showPlot , 

40 hold off, plot([0 5*alphal] ,[f0(l) f 0(l)+5*alphal*g0] , 'b- ' ) , hold on 

41 plot(0,f0(l) ,'bo'), 

42 s=-alphal :alphal/10:5*alphal ; 

43 end 

44 

45 while n<9, 

46 '/, Evaluate function at x+alphal . 

47 x = xO+alphal ; 

48 fl = eval (fcall); n = n+1 ; 

49 if showSteps, disp(['l: ' ,num2str (alphal) , ' ' ,num2str(f 1) , ' ',how]), end 

50 if showPlot, plot (alphal, f 1(1) , 'bo') , end 

51 

52 '/, Fit quadratic function to the points x0,xl. f = a*x~2 + b*x + c. 

53 c = f0(l); b = gO; a = (f 1 (l)-b*alphal-c)/(alphal*alphal) ; 

54 

55 if showPlot, plot (s ,polyval( [a,b,c] ,s) , 'r- ' ) , end '/, Fitted curve 

56 

57 '/, Determine next jump based on the quadratic fit. 

58 if (a>eps) , '/, Jump forward but not too much. 

59 alphaq = -b/(2*a) ; 

60 if (alphaq<5*alphal) & (alphaq~=alphal) , 

61 alpha2 = alphaq; 

62 how = [how , 'quadratic fit; '] ; 

63 else 

64 alpha2 = 5*alphal ; 

65 how = [how , 'limited jump; '] ; 

66 end 

67 

68 elseif (a<eps) , '/, Optimum is maximum, jump forward. 

69 alpha2 = 5*alphal ; 

70 how = [how ,' Increase step size (q) ; '] ; 
71 

72 end 
73 

74 '/, Evaluate function at x+alpha2 . 

75 x = xO + alpha2 ; 

76 f2 = eval (fcall); n = n+1; 

77 if showSteps, disp(['2: ' ,num2str (alpha2) , ' ' ,num2str(f 2(1) ) , ' '.how]), end 

78 if showPlot, plot (alpha2,f 2(1) , 'x' ) , end '/, Chosen point 
79 

80 '/, Check for adequate solution 

81 '/. if (f2(l)<min(f0(l),fl(l))) & (alpha2<alphal) , 

82 '/, alpha = alpha2 ; 

83 '/. f = f2; 

84 '/, break '/, Normal exit 

85 '/. end 

86 

87 '/, Try to fit a cubic function, f = a*x~3 + b*x~2 + c*x + d. 

88 d = f0(l) ; c = gO; 

89 alphamat = [alphal . "3 alphal . ~2 ; alpha2 . "3 alpha2.~2]; 

90 if rcond (alphamat )<eps , '/, Cubic solution is invalid, use quadratic 
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91 alpha = alpha2 ; 

92 f = f2; 

93 return '/, Normal exit 

94 end 

95 ab = alphamat\[fl(l)-c*alphal-d;f2(l)-c*alpha2-d] ; 

96 a = ab(l); b = ab(2) ; 

97 

98 if showPlot , plot (s ,polyval( [a,b,c ,d] ,s) , 'g- ' ) , end '/, Fitted curve 

99 

100 '/, Check descriminent 

101 del = b*b-3*a*c; 

102 if (del>0) & abs(a)>eps, '/, Solve for roots. 

103 alpha3 = (-b+sqrt (del) )/(3*a) ; 

104 elseif abs(a)<eps '/, Function looks VERY quadratic. 

105 alpha3 = alpha2 ; 

106 else 

107 alpha3 = inf ; '/, Flag invalid fit. 

108 end 

109 

110 if alpha3<0 , '/, Reduce step size and start over. 

111 alphal = min(alphal ,alpha2)/2 ; 

112 how = [how, 'Reduce step size (c); '] ; 

113 if showSteps, disp(['3: ' ,num2str(alpha3) , ' ',how]), end 

114 

115 elseif alpha3==inf , '/, Increase step size and start over. 

116 alphal = 5*max (alphal ,alpha2) ; 

117 how = [how ,' Invalid cubic fit; '] ; 

118 if showSteps, disp(['3: ' ,num2str(alpha3) , ' ',how]), end 

119 

120 else 

121 if (alpha3-alpha2) /alphal < .01, '/, Extra function evaluation not necessary 

122 f3 = inf; 

123 else 

124 how = [how, 'Cubic fit; '] ; 

125 '/, Evaluate function at x+alpha3. 

126 x = xO + alpha3; 

127 f3 = eval(fcall); n = n+1 ; 

128 if showSteps, disp(['3: ' ,num2str(alpha3) , ' ' ,num2str (f3(l) ) , ' '.how]), end 

129 if showPlot, plot (alpha3 ,f3(l ) , 'x ' ) , end '/, Chosen point 

130 end 

131 

132 '/, Check if minimum is bracketed. 

133 ftest = [fl(l),f2(l),f3(D] ; 

134 i = f ind(min(ftest )==ftest ) ; 

135 alpha = eval( ['alpha' ,int2str(i(l) )]) ; 

136 f = eval(['f ' ,int2str(i(l))]); 

137 

138 if (f(l)<fO(D) & (abs(alpha)>tol), 

139 return; '/, Normal exit 

140 

141 elseif abs(alpha)<tol , '/, GRAD may not be accurate. 

142 x = xO + tol; 

143 f = eval(fcall); n = n+1; 

144 alpha = 10*rand(l)*tol ; 

145 how = [how, 'Random Jump; '] ; 

146 return 

147 

148 else 

149 '/.if alpha==min( [alphal ,alpha2 ,alpha3] ) , 

150 alphal = alpha/2; 

151 how = [how, 'Reduce step size (b) ; '] ; 

152 

153 '/. else 

154 '/, alphal = alpha; 

155 '/, how = [how, 'Choose smallest; '] ; 

156 end 

157 

158 end 

159 

160 end '/, while 

161 

162 '/, Abnormal exit . No improvement in cost . Use min value for alpha 

163 f = f0; 

164 alpha = 0; 

165 disp( 'WARNING : Exceeded 9 line search attempts.') 

B.2.6 filter2d.m 

1 function x = filter2d(a, stencil , resize) 

2 '/.FILTER2D Two dimensional computational stencil filtering. 

3 '/. 

4 '/. X = FILTER2DC A, STENCIL) returns X which is the result of 

5 '/, applying the computational molecule STENCIL to the matrix A. 
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6 '/. 

7 '/. If SIZE(A) is ma-by-na and SIZE (STENCIL) is ms-by-ns then 

8 '/. SIZE(X) is (ma+ms-l)-by-(na+ns-l) . 

9 '/. 

10 '/. FILTER2D (A, STENCIL) is the same as C0NV2 (A.ROT90 (STENCIL ,2) ) . 

11 

12 '/. Clay M. Thompson 1-15-91 

13 

14 [ms,ns] = size(stencil) ; 

15 [ma,na] = size(a); 

16 

17 '/, This calculation is the same as conv2(a,rot90(stencil ,2) ) 

18 '/, but is faster. 

19 if 0, 

20 x = zeros (ms+ma-1 ,ns+na-l) ; 

21 for i=l:ms 

22 for j=l :ns 

23 w = stencil(ms-i+l ,ns-j+l) ; 

24 if w~=0, 

25 x = x + [zeros(i-l , j-1) ,zeros(i-l ,na) ,zeros(i-l ,ns-j ) ; 

26 zeros (ma, j-1) ,w*a,zeros(ma,ns-j ) ; 

27 zeros (ms-i ,j-l) , zeros (ms-i ,na) , zeros (ms-i ,ns-j )] ; 

28 end 

29 end 

30 end 

31 

32 else 

33 x = conv2(a,rot90(stencil ,2)) ; 

34 end 

35 

36 if nargin==3, '/, Return the central (valid) part. 

37 rows = ms-1 + [l:ma-ms+l]; 

38 cols = ns-1 + [l:na-ns+l]; 

39 x = x(rows ,cols) ; 

40 end 

41 



B.2.7 cfilter2d.m 

1 function c=cf ilter2d(a,b ,grad) 



6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 



,CFILTER2D Filter bicubic approximation of array. 

C = CFILTER2D(A,B) applies the filter B to the array A where 
A is interpolated using bicubic interpolation. 

C = CFILTER2D(A,B, 'grad') returns the convolution of B with A 
for a gradient calculation. 

The cubic approximation is used extrapolate a 1 pixel border 
around the array A. 

See also: FILTER2D. 

Clay M. Thompson 4-30-92 

ma,na] = size(a); '/, Size of array. 
mb,nb] = size(b); '/, Filter size. 



20 if nargin==2, 

21 if nb>l, 

22 a = [3*a(: ,l)-3*a(: ,2)+a(: ,3) ,a,3*a(: ,na)-3*a(: ,na-l)+a(: ,na-2)] ; 

23 end 

24 if mb>l, 

25 a = [3*a(l, : )-3*a(2, : )+a(3, : ) ;a;3*a(ma, : )-3*a(ma-l , : )+a(ma-2, : )] ; 

26 end 

27 c = f ilter2d(a,b, 'resize ') ; 

28 else 

29 c = f ilter2d(a,b) ; 

30 if nb>l, 

31 c(:,2:4) = c(:,2:4) + c(:,l)*[3 -3 1]; 

32 c( : ,na+nb-4:na+nb-2) = c( : ,na+nb-4 :na+nb-2) + c( : ,na+nb-l)* [1 -3 3]; 

33 c = c( : ,2 :nb+na-2) ; 

34 end 

35 if mb>l , 

36 c(2:4,:) = c(2:4,:) + [3 ;-3 ; l]*c(l , : ) ; 

37 c(ma+mb-4:ma+mb-2, : ) = c(ma+mb-4 :ma+mb-2 , : ) + [1 ; -3 ;3]*c (ma+mb-1 , : ) ; 

38 c = c(2 :mb+ma-2 , : ) ; 

39 end 

40 end 
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B.2.8 hbasis.m 



12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
38 
39 
40 
41 
42 
43 
44 
45 
46 
47 
48 
49 
50 
51 
52 
53 
54 
55 
56 
57 
58 
59 
60 
61 
62 
63 
64 
65 
66 
67 
68 
69 
70 
71 
72 
73 
74 
75 
76 
77 
78 
79 
80 
81 



function y = hbasis(x ,n,oper) 
HBASIS Map from hierarchical basis representation to nodal 
representation. 



Y = HBASIS(X,N) or Y = HBASIS 
hierarchical basis with N lev 

Y = HBASIS(X,N, 'trans') maps 
hierarchical basis with N lev 

Y = HBASIS(X,N, 'inv') maps fr 
hierarchical basis with N lev 

Y = HBASIS(X,N, 'tinv') maps f 
hierarchical basis with N lev 



Reference : 
Functions" 



"Fast Surface Inte 
Richard Szeliski 



(X , N, 'nodal ' ) maps X from a 
els to the nodal basis. 

from the nodal basis to the 
els using the adjoint map. 

om the nodal basis to the 
els . 

rom the nodal basis to the 
els using the adjoint map. 

rpolation Using Heirarchical Basis 
IEEE PAMI, Vol 12, No. 6, June 1990. 



Clay M. Thompson 4-2-91 
Revised 6-10-91 by CMT 

error(nargchk(2 ,3 ,nargin) ) ; 

if narginO , code = 'no'; else co 



de = oper(l:2); end 

x is compatible with n levels. 



'/, Check to make sure the size of 

[mm,nn] = size(x); 

incr = 2~(n-l) ; 

if mm<incr | nn<incr, 

error(['Each dimension of X must be larger than ' , int2str(incr) , 

' for ' ,int2str(n) , ' levels.']); 
end 

if code=='no', '/, oper== ' nodal ' 
y = zeros(mm,nn) ; 

rows = l:incr:mm; cols = l:incr:nn; 
y(rows,cols) = x(rows ,cols) ; 
for level=(n-l) :-l : 1 , 

xrows = 2*rem(mm-l ,incr) >= i 

incr = incr/2 ; 

row2 = l:incr:mm; col2 = 1 

y(row2,col2) = hb(y (rows , cols 

y(rows,cols) = y(rows,cols) - 

rows = row2 ; cols = col2 ; 
end 



ncr; xcols = 2*rem(nn-l ,incr) >= incr; 

ncr :nn ; 

), xrows , xcols) + x(row2 ,col2) ; 
x(rows ,cols) ; '/, Remove extra term 



elseif code=='tr', '/, oper=='trans 

y = zeros(mm,nn) ; 

incr = 1 ; 

rows = l:incr:mm; cols = 1 : inc 

y(rows,cols) = x(rows ,cols) ; 

for level=l : (n-1) , 
incr = 2*incr; 
row2 = l:incr:mm; col2 = 1 
y(row2,col2) = y(row2,col2) + 
rows = row2 ; cols = col2 ; 

end 

elseif code=='in', '/, oper=='inv', 

incr = 1 ; 

y = zeros(mm,nn) ; 

rows = l:incr:mm; cols = 1 : inc 

y(rows,cols) = x(rows ,cols) ; 

for level=l : (n-1) , 
incr = incr*2 ; 
row2 = l:incr:mm; col2 = 1 
xrows = 2*rem(mm-l ,incr) >= 
y(rows,cols) = y(rows,cols) 
y(row2,col2) = y(row2,col2) + 
rows = row2 ; cols = col2 ; 

end 



ncr :nn ; 
hbt (y(rows,cols)) ; 



ncr :nn ; 

ncr; xcols = 2*rem(nn-l ,incr) >= incr; 

hb(y(row2,col2) , xrows , xcols) ; 

x(row2 ,col2) ; '/, Add back extra term 



elseif code=='ti', '/, oper== 'tinv' , 

y = zeros(mm,nn) ; 

rows = l:incr:mm; cols = l:incr:nn 

y = x; 

for level=(n-l) :-l : 1 , 
incr = incr/2 ; 
row2 = l:incr:mm; col2 = 1 
y(rows,cols) = y(rows,cols) 



ncr :nn ; 
hbt(y(row2,col2)) ; 
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82 rows = row2 ; cols = col2 ; 

83 end 

84 

85 else 

86 error(['The operation ',oper,' is invalid.']); 

87 end 

88 
89 



B.2.9 hb.m 



1 function y=hb(x, xtrarows ,xtracols) 

2 '/,HB Interpolate a 2-D function on a rectangular grid using 

3 '/, Hierarchical Basis functions. 

4 '/. 

5 '/, Y = HB(X) returns X interpolated to the next level. If X is 

6 */. m-by-n then Y is (2*m-l)-by-(2*n-l ) . 

7 

8 */. Clay M. Thompson 4-2-91 

9 

10 colmask = [.5 .5]; 

11 rowmask = [.5;. 5]; 

12 

13 if nargin<2 , xtrarows = 0; end 

14 if narginO , xtracols = 0; end 

15 

16 [m,n] = size(x) ; 

17 if xtrarows | xtracols, 

18 x = [x , zeros (m, xtracols) ; zeros (xtrarows , n) , zeros (xtrarows , xtracols)] ; 

19 end 

20 

21 [m,n] = size(x) ; 

22 mm = 2*m-l ; nn = 2*n-l ; 

23 oddrows = l:2:mm; 

24 evenrows = 2:2:mm-l; 

25 oddcols = l:2:nn; 

26 evencols = 2:2:nn-l; 

27 

28 y = zeros(mm,nn) ; 

29 y (oddrows , oddcols) = x; 

30 y (oddrows , evencols) = f ilter2d(x, colmask, 'resize ') ; 

31 y (evenrows ,: ) = filter2d(y (oddrows ,:) , rowmask, 'resize ') ; 

32 

33 if xtrarows | xtracols, 

34 y = y (1 :mm-xtrarows , 1 :nn-xtracols) ; 

35 end 

B.2.10 hbt.m 

1 function y = hbt(x) 

2 '/.HBT Decimate a 2-D function on a rectangular grid using Hierarchical Basis 

3 '/, functions. 

4 '/. 

5 '/, Y = HBT(X) returns X decimated to the next coarsest level. 

6 '/. If X is m-by-n then Y is ((m+l)/2)-by-( (n+l)/2) . 

r '/. 

8 '/, Note (m+1) and (n+1) must be divisible by 2. 

9 

10 '/. Clay M. Thompson 4-2-91 
11 

12 [m,n] = size(x) ; 

13 rows = [l:2:m]; cols = [l:2:n]; 

14 a = [0 ,zeros(l ,n) ,0 ;zeros(m, 1) ,x, zeros (m, 1) ;0,zeros(l ,n) ,0] ; 

15 

16 y = a(rows, cols) + 2*a(rows, cols+1) + a(rows, cols+2) + ... 

17 2*a(rows+l ,cols) + 2*a(rows+l , cols+2) + ... 

18 a(rows+2 ,cols) + 2*a(rows+2 , cols+1) + a(rows+2 , cols+2) ; 

19 y = y/4; 

20 
21 



B.2.11 interpx.m 

l function F=interpx(E,x) 



'/.INTERPX Linear interpolation in the x-direction. 

•/. 

'/, F = INTERPX(E,X) returns a matrix F containing the values of E 

'/, at the points X. The matrix X must have the same number of rows as 

'/, E. The values in the matrix X must be between 1 and N, where N is 

'/, the number of columns of E. The value NaN will be returned where 

'/, this is not the case. 



154 



APPENDIX B. M-FILE LISTINGS 



10 */. Clay M. Thompson 10-17-90 



[m,n] = size(E) ; 

[mx,nx] = size(x); 

if mx~=m, error('X must have the same number of rows as E.'); end 

'/, Compute nearest x position, xlow. 

xlow = f loor (x) ; 

out = (x==n) ; 

nout = sum(out(:)); 

if nout>0, xlow(out) = (n-l)*ones(nout , 1) ; end 

'/. Check for out of range values of x and set to 1 

out = (x<l) | (x>n) ; 

nout = sum(out(:)); 

if nout>0, x(out) = ones (nout , 1) ; xlow(out) = ones(nout , 1) ; end 

'/, Determine index into matrix elements 

'/, Note: y=[l :m] '*mones(l :nx) ; 

elem = [1 :m] '*mones(l :nx) + (xlow-l)*m; 

F = E(elem) + (x-xlow) .* (E(elem+m)-E(elem) ) ; 

'/, Set values of F where x is out of range to NaN. 
if any(out(:))> F(out) = NaN*ones(nout , 1) ; end 



B.2.12 domain2d.m 



i 

2 

3 

4 

5 
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8 
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22 

23 



function [x,y] = domain2d(x ,y) 
D0MAIN2D Generate X and Y arrays for 3-d plots. 

[XX, YY] = D0MAIN2D(X,Y) transforms the domain specified by vectors 
X and Y into arrays XX and YY that can be used for the evaluation 
of functions of two variables with 3-d mesh plots. For example, 
to evaluate the function x*exp(-x~2-y~2) over the range -2 < x < 2 
-2 < y < 2, 



[x,y] = meshdom(-2: .2:2, -2: 
z = x .* exp(-x.~2 - y.~2); 
mesh(z) 



.2:2) ; 



'/. 



J.N. Little 12-2-85 
Revised 20-May-90 , LS . 
Copyright (c) 1985, 1986, 



1990 by the MathWorks, Inc. 



nx = length(x) ; 

ny = length(y) ; 

x = x( : ) . ' ; '/, make sure 

x = x(ones(ny, 1),:); 

y = y(ny:-l:l) ; y = y(:); 

y = y(:,ones(l, nx)); 



x is a row vector 



'/, make sure y is a column vector 



B.2.13 icubic.m 
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function F=icubic(x ,y ,u) 
ICUBIC Cubic Interpolation of a 1-D function. 

F=ICUBIC(Y,X1) returns the value of the 1-D function Y at the 
points XI using cubic interpolation. length(F)=length(Xl ) . XI is 
an index into the vector Y. Y is the value of the function 
evaluated uniformly on a interval. If Y is a matrix, then 
the interpolation is performed for each column of Y. 

If Y is of length N then XI must contain values between 1 and N. 
The value NaN is returned if 
this is not the case . 

F = ICUBIC(X , Y,X1 ) uses the vector X to specify the coordinates 
of the underlying interval. X must be equally spaced and 
monotonic . XI must lie within the coordinates in X. 

See also: ILINEAR. 

Clay M. Thompson 7-4-91 

Based on "Cubic Convolution Interpolation for Digital Image 
Processing", Robert G. Keys, IEEE Trans, on Acoustics, Speech, and 
Signal Processing, Vol. 29, No. 6, Dec. 1981, pp. 1153-1160. 

if nargin==2 , '/, No X specified. 
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27 u = y ; y = x ; 

28 '/, Check for vector problem. If so, make everything a column vector. 

29 if min(size(y) )==1 , y = y(:); end 

30 if min(size(u) )==1 , u = u(:); end 

31 [nrows ,ncols] = size(y); 

32 

33 elseif nargin==3 , '/, X specified. 

34 '/, Check for vector problem. If so, make everything a column vector. 

35 if min(size(y) )==1 , y = y(:); end 

36 if min(size(x) )==1 , x = x(:); end 

37 if min(size(u) )==1 , u = u(:); end 

38 [nrows ,ncols] = size(y); 

39 '/, Scale and shift u to be indices into Y. 

40 if (min(size(x) )~=1 ) , error('X must be a vector.'); end 

41 x = x( : ) ; 

42 [m,n] = size(x); 

43 if m ~= nrows, 

44 error ('The length of X must match the number of rows of Y.'); 

45 end 

46 u = 1 + (u-x(l))*((nrows-l)/(x(m)-x(l))) ; 

47 

48 else 

49 error ( 'Wrong number of input arguments.'); 

50 end 

51 

52 if nrows<3, error('Y must have at least 3 rows.'); end 

53 [m,n] = size(u); 

54 if n==l , u = u*ones(l ,ncols) ; [m,n] = size(u); end '/, Expand u 

55 if n~=ncols , error('The number of columns in XI and Y must match.'); end 

56 

57 '/, Check for out of range values of u and set to 1 

58 uout = f ind((u<l) | (u>nrows)) ; 

59 nuout = length(uout ) ; 

60 if nuout>0, u(uout) = ones(nuout , 1) ; end 

61 

62 '/, Interpolation parameters 

63 s = (u - floor(u)); 

64 u = f loor(u) ; 

65 d = find(u==nrows) ; if length(d)>0, u(d) = u(d)-l; s(d) = s(d)+l; end 

66 

67 '/, Expand y so interpolation is valid at the boundary. 

68 y = [3*y(l , : )-3*y(2, : )+y(3, : ) ;y ;3*y (nrows , : )-3*y (nrows- 1 , : )+y(nrows-2, : )] ; 

69 nrows = nrows + 2 ; 

70 

71 '/, Now interpolate using computationally efficient algorithm. 

72 s2 = s.*s; s3 = s.*s2; 

73 ndx = u+ones(m,l)+ [0:n-l]*nrows ; 

74 F = y(ndx) .*(-s3+2*s2-s) + y (ndx+1 ) .*(3*s3-5*s2+2) + ... 

75 y(ndx+2) . *(-3*s3+4*s2+s) + y (ndx+3) .*(s3-s2) ; 

76 F = F/2; 

77 

78 '/, Now set out of range values to NaN . 

79 if nuout>0, F(uout) = NaN*ones (nuout ,1 ) ; end 

80 



B.2.14 dcubicx.m 

1 function G=dcubicx(x,y ,u) 
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'/.DCUBICX Derivative of Cubic Interpolation of a 1-D function w.r.t. x.. 

G = DCUBICX(Y,X1) returns the derivative of the 1-D function Y=f(x) 
at the points XI using cubic interpolation. length(G)=length(Xl) . 
XI specifies the points originally used for the interpolation. 

If Y is of length N then XI must contain values between 1 and N. 
The value NaN is returned if this is not the case . 

G = DCUBICX(X ,Y ,X1) uses the vector X to specify the coordinates 
for Y as for ICUBIC. 

See also: ICUBIC. 

Clay M. Thompson 7-18-91 

Based on "Cubic Convolution Interpolation for Digital Image 
Processing", Robert G. Keys, IEEE Trans, on Acoustics, Speech, and 
Signal Processing, Vol. 29, No. 6, Dec. 1981, pp. 1153-1160. 



20 
21 

22 if nargin==2 , '/, No X specified. 

23 u = y ; y = x; 

24 '/, Check for vector problem. If so, make everything a column vector 
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25 if min(size(y) )==1 , y = y(:); end 

26 if min(size(u) )==1 , u = u(:); end 

27 [nrows ,ncols] = size(y); 

28 

29 elseif nargin==3 , '/, X specified. 

30 '/, Check for vector problem. If so, make everything a column vector. 

31 if min(size(y) )==1 , y = y(:); end 

32 if min(size(x) )==1 , x = x(:); end 

33 if min(size(u) )==1 , u = u(:); end 

34 [nrows ,ncols] = size(y); 

35 '/, Scale and shift u to be indices into Y. 

36 if (min(size(x) )~=1 ) , error('X must be a vector.'); end 

37 x = x( : ) ; 

38 [m,n] = size(x); 

39 if m ~= nrows, 

40 error ('The length of X must match the number of rows of Y.'); 

41 end 

42 u = 1 + (u-x(l))*((nrows-l)/(x(m)-x(l))) ; 

43 

44 else 

45 error ( 'Wrong number of input arguments.'); 

46 end 

47 

48 if nrows<3, error('Y must have at least 3 rows.'); end 

49 [m,n] = size(u); 

50 if n==l , u = u*ones(l ,ncols) ; [m,n] = size(u); end '/, Expand u 

51 if n~=ncols , error ('The number of columns in XI and Y must match.'); end 

52 

53 '/, Check for out of range values of u and set to 1 

54 uout = f ind((u<l) | (u>nrows)) ; 

55 nuout = length(uout ) ; 

56 if nuout>0, u(uout) = ones(nuout , 1) ; end 

57 

58 '/, Interpolation parameters 

59 s = (u - floor(u)); 

60 u = f loor (u) ; 

61 d = find(u==nrows) ; if length(d)>0, u(d) = u(d)-l; s(d) = s(d)+l; end 

62 

63 '/, Expand y so interpolation is valid at the boundary. 

64 y = [3*y(l , : )-3*y(2, : )+y(3, : ) ;y ;3*y (nrows , : )-3*y (nrows- 1 , : )+y(nrows-2, : )] ; 

65 nrows = nrows + 2 ; 

66 

67 '/, Now interpolate using computationally efficient algorithm. 

68 s2 = s.*s; 

69 ndx = u+ones(m,l)+ [0:n-l]*nrows ; 

70 G = y(ndx) .*(-3*s2+4*s-l) + y (ndx+1) . *(9*s2-10*s) + ... 

71 y(ndx+2) . *(-9*s2+8*s+l) + y (ndx+3) . *(3*s2-2*s) ; 

72 G = G/2; 

73 

74 '/, Now set out of range values to NaN . 

75 if nuout>0, G(uout) = NaN*ones (nuout ,1 ) ; end 

76 



B.2.15 dcubicz.m 

l function g=dcubicy (x,y ,u ,dJdf ) 



'/.DCUBIC Derivative of 1-D cubic interpolation w.r.t. Y. 

G = DCUBICY(Y,Xl,dJdF) computes the derivative of the cost 
function J(F) with respect to the underlying variables Y. The 
matrix dJdF is the derivative of the cost function with respect to 
the interpolated value F = ICUBIC(Y ,X1) . XI specifies the points 
originally used for the interpolation. 



See also: ICUBIC. 



10 7. G = DCUBICY(X,Y,X1 ,dJdF) uses the vector X to specify the 
coordinates for Y as for ICUBIC. 

12 
13 
14 

15 */. Clay M. Thompson 7-12-91 

16 

17 if nargin==3, '/, No X specified. 

18 dJdf =u; u=y; y=x; 

19 '/, Check for vector problem. If so, make everything a column vector. 

20 if min(size(y) )==1 , y = y(:); end 

21 if min(size(u) )==1 , u = u(:); end 

22 [nrows ,ncols] = size(y); 

23 

24 elseif nargin==4 , '/, X specified. 

25 '/, Check for vector problem. If so, make everything a column vector. 

26 if min(size(y) )==1 , y = y(:); end 
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27 
28 
29 
30 
31 
32 
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50 
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55 
56 
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100 
101 
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if min(size(x) )==1 , x = x(:); end 

if min(size(u) )==1 , u = u(:); end 

[nrows ,ncols] = size(y); 

'/, Scale and shift u to be indices into Y. 

if (min(size(x) )~=1 ) , error('X must be a vector.'); end 

x = x(:); 

[m,n] = size(x) ; 

if m ~= nrows, 

error('The length of X must match the number of rows of Y.'); 
end 
u = 1 + (u-x(l))*((nrows-l)/(x(m)-x(l))) ; 

else 

error('Wrong number of input arguments.'); 
end 

if nrows<3, error('Y must have at least 3 rows.'); end 

[m,n] = size(u) ; 

if n==l , u = u*ones(l ,ncols) ; [m,n] = size(u); end '/, Expand u 

if any (size(u) ~=size(dJdf ) ) , error('dJdF and XI must be the same size.'); end 

if n~=ncols , error('The number of columns in XI and Y must match.'); end 

if m<2 , '/, Expand u (so sums work) with an out of range value. 

u = [u ;zeros(l ,n)] ; [m,n] = size(u); 

dJdf = [dJdf ;zeros(l,n)] ; 
end 

'/, Check for out of range values of u and set to 1 

uout = f ind( (u<l) I (u>nrows) ) ; 

nuout = length(uout ) ; 

if nuout>0, u(uout) = ones(nuout , 1) ; end 

'/, Interpolation parameters 

s = (u - f loor(u) ) ; 

u = f loor(u) ; 

d = f ind(u==nrows) ; if length(d)>0, u(d) = u(d)-l; s(d) = s(d)+l; end 

'/, Compute terms for gradient 
s2 = s.*s; s3 = s.*s2; 
tO = dJdf .*(-s3+2*s2-s) ; 
tl = dJdf .*(3*s3-5*s2+2); 
t2 = dJdf .*(-3*s3+4*s2+s) ; 
t3 = dJdf .*(s3-s2); 
clear s s2 s3 

'/, Set out of range terms to zero, 
if nuout >0 , 

tO(uout) = zeros(nuout , 1) ; 

tl(uout) = zeros(nuout , 1) ; 

t2(uout) = zeros (nuout , 1) ; 

t3(uout) = zeros(nuout , 1) ; 

dJdf(uout) = zeros(nuout , 1) ; 
end 

'/, Form sums for each z value . 
tOsum = zeros(nrows ,ncols) ; 
tlsum = zeros(nrows ,ncols) ; 
t2sum = zeros(nrows ,ncols) ; 
t3sum = zeros(nrows ,ncols) ; 
for k=l: nrows, 

elem = find(u==k); 

temp = zeros(m,n); 

temp(elem) = tO(elem) ; t0sum(k,:) = sum(temp) 

temp(elem) = tl(elem); tlsum(k,:) = sum(temp) 

temp(elem) = t2(elem); t2sum(k,:) = sum(temp) 

temp(elem) = t3(elem) ; t3sum(k,:) = sum(temp) 
end 
clear tO tl t2 t3 



'/, Add terms from boundary conditions 



tlsum(l , 
t2sum(l , 
t3sum(l , 

t2sum(m-l , 
tlsum(m-l , 
t0sum(m-l , 



tlsum(l , 
= t2sum(l, 
= t3sum(l, 



) + 3*t0sum(l , : 



3*t0sum(l 
t0sum(l , : ) 



) = t2sum(m-l , 
) = tlsum(m-l , 
) = t0sum(m-l , 



) + 3*t3sum(m-l , 
) - 3*t3sum(m-l , 
) + t3sum(m-l , : ) 



'/, Now combine to form gradient . 

g = zeros(nrows ,ncols) ; 

g(l,:) = t0sum(2,:) + tlsumd,:); 

g(2,:) = t0sum(3,:) + tlsum(2,:) + t2sum(l , : ) ; 

for k=3:nrows-l, 

g(k,:) = t0sum(k+l,:) + tlsum(k,:) + t2sum(k-l,:) + t3sum(k-2 , : ) ; 
end 
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111 g(nrows,:) = tlsum(nrows , : ) + t2sum(nrows-l , : ) + t3sum(nrows-2 , : ) ; 

112 

113 g = g/2; 

114 
115 
116 



Appendix C 

Test Surface Descriptions 



This appendix contains information on how to create the crater and hill test surfaces 
and associated test images. 

C.l Crater Surface 

The crater on the flat plane is formed by intersecting two spheres, one with radius 
fO and the other with radius 9.48, with a flat plane. The routines makepair4 and 
crater_depth were used to create the crater images. 

C.l.l makepair4.m 

1 '/, Make stereo pair (this is the best routine) 

2 clear 

3 DEPTH = 'crater_depth' 

4 '/.DEPTH = 'sphere_depth' 

5 '/.GRAD = 'crater_grad' 

6 n = 65 ; '/, The size of image . 

7 mu = 1 ; '/, Optimization parameter 

8 

9 b = 500; 
10 f = -1000; '/, We are working in the same coordinates as the surface for x & y. 
11 

12 '/, Camera coordinate calibration 

13 zO = eval( [DEPTH ,'(0,0)']) 

14 vO = 0; 

15 vl = vO + b*f /zO/2 

16 v2 = vO - b*f /zO/2 
17 

18 gamma = f *b/2 ; 

19 delta = zO/f ; 

20 if 1 '/, Case a {Hard case} 

21 psl = .1; qsl = .1; 

22 ps2 = -. 1 ; qs2 = .1 ; 

23 prefix = ['cr' ,int2str(n) , 'a'] ; 

24 else '/, Case b {Easy case} 

25 psl = .2; qsl = -.5; 

26 ps2 = -.3; qs2 = .1 ; 

27 prefix = ['cr' ,int2str(n) , 'b'] ; 

28 end 

29 '/.psl = .5; qsl = -.6; 

30 '/.ps2 = .5; qs2 = -.6; 

31 '/.psl = 0; qsl = 0; 

32 

33 [x,y] = domain2d(-12:24/(n+l):12,-12:24/(n+l):12 ); 

34 

35 '/, Image parameters 
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36 deltax = x(l,2)-x(l ,1) ; 

37 params = [f/deltax '/, Focal length for unity spacing grid 

38 b '/, baseline distance 

39 zO '/, Calibration depth 

40 psl ;qsl ;ps2 ;qs2 '/, Light positions 

41 (x(2 ,2)+vl)/deltax-l '/, x-calibration as offset from E(l,l). 

42 y (2 ,2)/deltax-n '/, y-calibration as offset from E(l,l). 

43 (x(2 ,2)+v2)/deltax-l '/, x-calibration as offset from E(l,l). 

44 y (2 ,2)/deltax-n '/, y-calibration as offset from E(l,l). 

45 1 '/, Grid spacing 

46 ] ; 

47 del = zO*(max(x(: ))-min(x(: )))/f/(n+l) ; 

48 

49 z = eval([DEPTH, >(x,y)>]); 

50 eqmesh(z ,del) 

51 titleOTrue Surface') 

52 

53 '/, Determine depth function for left image (1) 

54 

55 err = 1; m = 0; 

56 zl = zO*mones(x); 

57 gamma/zO 

58 while (err> 1 .e-5)ft(m<25) , 

59 ztemp = eval([DEPTH, ' (x+vl-gamma. /zl-vO ,y) '] ) ; 

60 err = sum( (ztemp( : )-zl ( : ) ) . ~2 ) 

61 zl = zl + mu*(ztemp-zl) ; 

62 m = m + 1 ; 

63 end 

64 eqmesh(zl ,del) 

65 title ('Left Surface') 

66 

67 '/, determine depth function for right image (2) 

68 err = 1; m = 0; 

69 z2 = z0*mones(x); 

70 while (err>l.e-5)&(m<25) , 

71 ztemp = eval([DEPTH, ' (x+v2+gamma./z2-v0,y) '] ) ; 

72 err = sum( (ztemp( : )-z2( : ) ) . ~2 ) 

73 z2 = z2 + mu*(ztemp-z2) ; 

74 m = m + 1 ; 

75 end 

76 eqmesh(z2 ,del) 

77 title ('Right Surface') 
78 

79 hx = [-1 l]/2; 

80 hy = [l;0;-l]/2; 

81 rows = 2:n+l; 

82 cols = 2:n+l; 

83 

84 '/, Right Image 

85 zx = filter2d(zl ,hx , 'resize ' )/deltax ; 

86 zy = filter2d(zl ,hy , 'resize ' )/deltax ; 

87 pi = f*zx(rows ,:)./( (x(rows ,cols)+vl) . *zx(rows ,: )+zl (rows , cols) ) ; 

88 ql = f*zy (: ,cols) . /(y (rows ,cols) . *zy( : ,cols)+zl (rows ,cols) ) ; 

89 El = rmap(pl ,ql ,psl ,qsl) ; 

90 

91 '/, Left Image 

92 zx = filter2d(z2 ,hx, 'resize ' )/deltax; 

93 zy = filter2d(z2 ,hy, 'resize ' )/deltax; 

94 p2 = f*zx(rows ,:)./( (x(rows ,cols)+v2) . *zx(rows ,: )+z2(rows , cols) ) ; 

95 q2 = f *zy ( : ,cols) . /(y (rows ,cols) . *zy ( : ,cols)+z2(rows ,cols) ) ; 

96 E2 = rmap(p2,q2,ps2 ,qs2) ; 

97 

98 '/, Center Global coordinate system 

99 zx = filter2d(z ,hx, 'resize ' )/deltax ; 

100 zy = filter2d(z ,hy, 'resize ' )/deltax ; 

101 p = f *zx(rows ,:). /(x(rows , cols) .*zx(rows ,: )+z(rows , cols) ) ; 

102 q = f *zy( : ,cols) . /(y (rows , cols) .*zy( : ,cols)+z(rows , cols) ) ; 

103 E = rmap(p , q,ps2 ,qs2) ; 

104 

105 clg, 

106 subplot (121) , eqmesh(zl ,del) , 

107 subplot (122) , eqmesh(z2,del) 

108 subplot (111) , titleClmage contours') 

109 

110 mimage ([El .zeros (n,l) ,E2] , [0 1]) 

111 clear pi ql p2 q2 zx zy p q E gamma err m n ztemp 

112 clear psl ps2 qsl qs2 zO bO b f delta hx hy rows cols vl v2 

113 clear mu vO x y deltax 

114 
115 
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C.1.2 crater .depth, m 



16 
17 
18 
19 
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function [z ,rho]=depth(x ,y) 

*/,CRATER_DEPTH Depth function for crater on a flat plane. 

'/. 

'/. Z = CRATER_DEPTH(X,Y) . X and Y must be of the same size. 

'/. 

'/. Crater is defined between -10 <= x,y <= 10. 

zO = -1000; 

rl = 10; r2 = 9.48; 

dl = find(x.*x + y.*y <= 0.75*rl*rl); 

d2 = find(x.*x + y.*y <= 0.75*r2*r2); 

[m,n] = size(x) ; 

z = zO*ones(m,n) ; 

z(dl) = z(dl) + sqrt(rl*rl - x(dl).*x(dl) - y (dl) . *y (dl) )-5 ; 

z(d2) = min(z(d2), zO+12-sqrt (r2*r2 - x(d2).*x(d2) - y(d2) .*y(d2) )) ; 

'/.z = zO + (z-z0)/2; '/. Reduce height by two 

if nargout>l, 

*/. Albedo (rho) 

rho = mones (z) ; 

d = f ind(abs(x-y)<3) ; 

rho(d) = 0.7*mones(d) ; '/, Dark strip across diagonal 
end 



C.2 Hill Surface 

A fractal based method was used to create the underlying data matrix for the hill 
surface. Given this underlying data matrix the hill surface is formed using bicubic 
interpolation to define the surface points. The data for the hill is shown in the matrix 
(DATA_MATRIX) below, 



1000 + 



0.7210 


1.0934 


1.1456 


0.7672 


0.5582 


0.6868 


1.0628 


0.9350 


0.4301 


1.2280 


0.4993 


0.5555 


0.8574 


0.3552 


0.8864 


1.0791 


0.2720 


0.5742 


0.6650 


0.3307 


0.2967 


0.7894 


0.2520 


0.5136 


0.2594 


1.1026 


1.0038 


0.4796 


0.8690 


0.4103 


0.3271 


0.5276 


1.0706 


1.1861 


0.4380 


1.1673 


1.1548 


0.3879 


0.2889 


0.3674 


0.8726 


0.7950 


1.3493 


0.4233 


0.3691 


0.2568 


0.4342 


0.6595 


0.5902 



(c.i; 



The routines makepair_data and data_depth were used to create the hill images. 
The mountain images were formed using a similar method based on a 33-by-33 data 
matrix that is too large to be presented here. 



C.2.1 makepair data.m 

1 '/, Make stereo pair based on data 

2 clear '/, clear all data 

3 global DATA_MATRIX 

4 DEPTH = 'data_depth' 

5 n = 65 ; '/, The size of image . 

6 nn = n+2 ; '/, The size of z. 

7 mu = 1 ; '/, Optimization parameter 

8 

9 if 1 , '/. Wrinkled surface (hill) 

10 load data_matrix 

11 b = 500; 

12 psl = -1 ; qsl = 1 ; 

13 ps2 = .3; qs2 = .1; 

14 prefix = ['w' ,int2str(n) , 'c'] ; 
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15 else '/, Lake (mountain) 

16 load lake_data 

17 zO = -1000; 

18 DATA_MATRIX = DATA_MATRIX + zO; 

19 b = 100; 

20 psl = .5; qsl = .5; 

21 ps2 = -.5; qs2 = 0; 

22 pref ix = ['Ik' ,int2str(n) , 'b'] ; 

23 end 

24 [md,nd] = size(DATA_MATRIX) ; 

25 

26 f = -1000; '/, We are working in the same coordinates as the surface for x & y. 

27 

28 '/, Camera coordinate calibration 

29 zO = eval( [DEPTH ,'(0,0)']) 

30 vO = -b*f/z0/2; 

31 vl = vO + b*f /zO/2 

32 v2 = vO - b*f /zO/2 

33 

34 gamma = f *b/2 ; 

35 delta = zO/f ; 

36 

37 [x,y] = domain2d([l:(md-l)/(nn-l):md]-(md+l)/2,[l:(nd-l)/(nn-l):nd]-(nd+l)/2); 

38 x = x*.9; y = y*.9; '/. Center 90'/. of matrix 

39 

40 '/, Image parameters 

41 deltax = x(l,2)-x(l ,1) ; 

42 params = [f/deltax '/, Focal length for unity spacing grid 

43 b '/, baseline distance 

44 zO '/, Calibration depth 

45 psl ;qsl ;ps2 ;qs2 '/, Light positions 

46 (x(2 ,2)+vl)/deltax-l '/, x-calibration as offset from (1,1). 

47 y (2 ,2)/deltax-n '/, y-calibration as offset from (1,1). 

48 (x(2 ,2)+v2)/deltax-l '/, x-calibration as offset from (1,1). 

49 y (2 ,2)/deltax-n '/, y-calibration as offset from (1,1). 

50 1 '/, Grid spacing 

51 ]; 

52 del = zO*(max(x(: ))-min(x(: )))/f/(n+l) ; 

53 

54 z = eval([DEPTH, '(x,y)']); 

55 eqmesh(z ,del) 

56 '/.mesh(z) 

57 titleCTrue Surface') 

58 

59 '/, Determine depth function for left image (1) 

60 

61 err = 1; m = 0; 

62 zl = zO*mones(x); 

63 gamma/zO 

64 while (err> 1 .e-5)ft(m<25) , 

65 ztemp = eval([DEPTH, ' (x+vl-gamma. /zl-vO ,y) '] ) ; 

66 err = sum( (ztemp( : )-zl ( : ) ) . ~2 ) 

67 zl = zl + mu*(ztemp-zl) ; 

68 m = m + 1 ; 

69 end 

70 eqmesh(zl ,del) 

71 title ('Left Surface') 
72 

73 '/, determine depth function for right image (2) 

74 err = 1 ; m = ; 

75 z2 = zO*mones(x); 

76 while (err>l.e-5)ft(m<25) , 

77 ztemp = eval([DEPTH, ' (x+v2+gamma./z2-v0,y) '] ) ; 

78 err = sum( (ztemp( : )-z2( : ) ) . ~2 ) 

79 z2 = z2 + mu*(ztemp-z2) ; 

80 m = m + 1 ; 

81 end 

82 eqmesh(z2 ,del) 

83 title ('Right Surface') 

84 

85 hx = [-1 l]/2; 

86 hy = [l;0;-l]/2; 

87 rows = 2:n+l; 

88 cols = 2:n+l; 

89 

90 '/, Right Image 

91 zx = filter2d(zl ,hx , 'resize ' )/deltax ; 

92 zy = filter2d(zl ,hy, 'resize ' )/deltax; 

93 pi = f*zx(rows ,:)./( (x(rows ,cols)+vl) . *zx(rows ,: )+zl (rows , cols) ) ; 

94 ql = f*zy (: ,cols) . /(y (rows ,cols) . *zy( : ,cols)+zl (rows ,cols) ) ; 

95 El = rmap(pl ,ql ,psl ,qsl) ; 
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97 '/, Left Image 

98 zx = filter2d(z2 ,hx , 'resize ' )/deltax ; 

99 zy = filter2d(z2 ,hy , 'resize ' )/deltax ; 

100 p2 = f *zx(rows ,:)./( (x(rows , cols)+v2) . *zx(rows , : )+z2(rows , cols) ) ; 

101 q2 = f *zy ( : ,cols) . /(y (rows ,cols) . *zy ( : ,cols)+z2(rows ,cols) ) ; 

102 E2 = rmap(p2,q2,ps2 ,qs2) ; 

103 

104 '/, Center Global coordinate system 

105 zx = filter2d(z ,hx, 'resize ' )/deltax ; 

106 zy = filter2d(z ,hy, 'resize ' )/deltax ; 

107 p = f *zx(rows ,:). /(x(rows , cols) .*zx(rows ,: )+z(rows , cols) ) ; 

108 q = f *zy( : ,cols) . /(y (rows , cols) .*zy( : ,cols)+z(rows , cols) ) ; 

109 E = rmap(p ,q,ps2 ,qs2) ; 
no 

ill clg,subplot(221) 

112 subplot (121) , eqmesh(zl ,del) 

113 subplot (122) , eqmesh(z2 } del) 

114 subplot (111) , title ( 'Image contours ' ) 

115 

116 mimage([El ,zeros(n,l) ,E2] , [0 1]) 

117 

118 clear pi ql p2 q2 zx zy p q E gamma err m n ztemp 

119 clear psl ps2 qsl qs2 zO bO b f delta hx hy rows cols vl v2 

120 clear mu vO x y md nd nn deltax 

121 

122 ztrue = z; 

123 



C.2.2 datadepth.m 



'/,DATA_DEPTH Depth function based on data matrix. 

•/. 

'/. Z = DATA_DEPTH(X,Y) . X and Y must be of the same size. 

'/. 

'/. [Z.RH0] = DATA_DEPTH(X,Y) returns albedo also. 

'/. 

'/. Relies on the global DATA_MATRIX to define surface. 

'/. 

'/. For DATA_MATRIX m-by-n, values in X must be between -(n-l)/2 and (n-l)/2, 

'/, values in Y must be between -(m-l)/2 and (m-l)/2. 



13 '/.global DATA_MATRIX 

14 [m,n] = size(DATA_MATRIX); 

15 [xx,yy] = domain2d(l :n,l :m) ; 

16 

17 x = x + (m+l)/2; y = y + (n+l)/2; 

18 

19 d = f ind(f loor(x)<l) ; x(d) = ones(length(d) , 1) ; 

20 d = find(x>n); x(d) = n*ones (length (d) , 1) ; 

21 d = f ind(floor(y)<l) ; y(d) = onesdength(d) , 1) ; 

22 d = f ind(y>m) ; y(d) = m*ones (length (d) , 1) ; 

23 

24 z = bicubic(xx,yy ,DATA_MATRIX,x,y) ; 

25 '/.z = blinear(xx,yy,DATA_MATRIX,x,y) ; 

26 

27 d = f ind(isnan(z) ) ; 

28 if length(d)>0, keyboard, end 

29 if length(d)>0, z(d) = onesdength(d) , l)*min(DATA_MATRIX( : ) ) ; end 

30 

31 if nargout>l, 

32 '/,d = f ind(abs(x-y)<l) ; 

33 zO = min(z( : )) ; 

34 d = find((z>.5+z0) & (z<.7+z0)); 

35 rho = mones(z); 

36 rho(d) = 0.7*mones(d) ; 

37 end 
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